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Turbulence Program for Propulsion Systems 

Tsan-Hsing Shih 


1. Introduction 

1.1 Background 

• The center for modeling of turbulence and transition (CMOTT) at NASA Lewis 
Research Center (LeRC) has been in existence for about four years. In the first 
three years, its main ax:tivities were developing and validating turbulence and 
combustion models for propulsion systems, in an effort to remove the deficiencies 
of existing models. Three workshops on computational turbulence modeling were 
held at LeRC (1991, 1993, 1994). 

• A peer review of the turbulence modeling activities at LeRC was held in Septem- 
ber, 1993. Seven peers from GE, P&W, RocketDyne, Cornell University, Uni- 
versity of California at Berkeley and NASA Ames conducted the review. The 
objective of the review was to assess the turbulence program at LeRC/CMOTT 
and to suggest the future direction of turbulence modeling activities for propul- 
sion systems. In September of 1994, a NASA-wide turbulence program peer 
review was held at NASA Ames Research Center (ARC) by sixteen peers from 
the aerospace industry, universities and other agencies. 

1.2 Recommendations from the 1993 peer review 

• “LeRC should spend substantial effort being responsive to industry’s current 
pressing perceived needs; this involves extensive discussion with industry during 
every phase of model development, analysis of industry’s problems, goal oriented 
model development, evaluation of models relative to industry’s intended applica- 
tion ...” 

• “LeRC has an obligation not only to respond to industry’s requests for help, but 
to play an autonomous, independent leadership role in providing models of the 
highest quality, . . ., which can be employed not only by the aircraft gas turbine 
and rocket industries but also by other industries . . .” 

• “In the present financial climate, industry does not have the resources to un- 
dertake model development and evaluation. LeRC’s help in this regard via the 
creation of its turbulence modeling effort, is, therefore, welcome from the indus- 
try’s standpoint.” 

• “It is important to work with the industry to evaluate the models and rank-order 
them by performance and cost in order to identify the most appropriate models 
for particular situations.” 

• There are many other useful suggestions and comments including collaboration 
with industry, joint programs, industry-wide workshops, etc. 
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1.3 Recommendations from the 1994 NAS A- wide turbulence peer review 

• NASA needs a right well defined turbulence program in order to help industry 
meet its needs. 

• The turbulence program needs a certain level of support, coordination, hnkage 
to industry needs, technology transfer, management and leadership. 

• The current progress is not reflective of real problems. 

Based on the peers’ comments and the current environment of NASA LeRC, we 
have planned and carried out a turbulence program for propulsion systems. The 
program is briefly explained in the following sections. 

2. Program goals at CMOTT 

• Develop reliable turbulence and bypass transition models and combustion models 
for complex flows in propulsion systems. 

• Integrate developed models into deliverable CFD tools for propulsion systems in 
collaboration with researchers at NASA LeRC and industry partners. 

3. Program approaches 

• Develop turbulence and combustion subprograms or modules for CFD users. 

• Develop collaboration program with NASA LeRC and industry partners to facil- 
itate technology transfer. 

• Conduct turbulence modeling research for propulsion systems: 

<> One-point moment closures for non-reacting flows, 

^ Scalar PDF method for turbulent reacting flows, 

0 Validation of existing and newly developed models. 

4. Development of turbulence and combustion subprogr ams (modules) 

4.1 Objective 

Build quick and efficient vehicles for turbulence and combustion technology trans- 
fer to CFD users in propulsion systems, for example, inlet/nozzle, turbomachinery 
and combustion, etc. 

4.2 The features of the turbulence module 

• It contains various turbulence models from which users can choose the model 
appropriate for the flows of interest. 

• It is self-contained, i.e., it contains its own solver for the turbulence model equa- 
tions. 

• It can be easily linked to industry CFD codes. 

4.3 The features of the PDF combustion module 

• It can be easily coupled with any existing industry flow codes. 

• It has a novel averaging scheme to reduce memory requirement. 

• It contains a general chemistry package. 

• There is a parallehzed workstation version. 
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4.4 NPARC Turbulence module 

NPARC is a CFD code widely used by NASA and industry. A turbulence sub- 
program has been developed for this code. 

• The models already built-in at the present time are: 

two mixing length models; the Chien k — s model; and the CMOTT k — s model. 

• The models to be built-in are: 

the CMOTT Reynolds stress algebraic equation model; a Reynolds stress trans- 
port equation model; and other models requested by CFD users. 

• A robust numerical solver is built-in for the model equations. 

4.5 V-STAGE turbomachinery turbulence module 

V-STAGE is an advanced turbomachinery code developed at NASA LeRC. 

• The built-in turbulence model is the CMOTT k — e model. 

5. Collaboration program and technology transfer 

5.1 Joint programs with NASA LeRC and industry 

• Preliminary programs with engine companies and others have been initiated (GE, 
P&W, RocketDyne, Naval Research Laboratories). 

• Joint programs with the oflBices of Advanced Subsonic Technology (AST), High 
Speed Research (HSR) and Numerical Propulsion System Simulation (NPSS) 
have been developed as ongoing programs with which industry partners will also 
be closely involved. 

5.2 Industry-wide workshop 

The industry-wide workshop will be a regular program held once every other year 
to: 

• release Lewis turbulence and combustion modules to CFD users; and 

• discuss the needs of industry and the state of the art in engineering turbulence 
modehng. 

6. Models developed at CMOTT 

Turbulence and combustion model development and validation are ongoing re- 
search activities. The following are a part of work done at CMOTT. 

• Isotropic eddy viscosity models: 

NASA TM 105633\ 1057682, 1062633, i06721^. 

• Reynolds stress and scalar flux algebraic equation models: 

NASA TM 1061163, i06513®, 106644^ IJNMF®. 

• Second moment transport equation models: 

NASA TM 105351^, 106469^®, 106681“, 105954“. 

• Multiple-scale models for compressible turbulent flows: 

NASA TM 106072“. 

• Bypass transition models: 
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NWTF^^. 

• PDF models for turbulent reacting flows: 
NASA TM IO66I415, AIAAI6, AIAA^^ 


PROGRAM SUMMARY 
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k - € Eddy Viscosity Turbulence Models: 
Development and Applications 

Z. Yang 


1. Motivation and Objective 

The purpose of this research is to develop k - e eddy viscosity models for the 
complex flows of engineering interest and to incorporate the resulting models into 
general purpose computational fluid dynamics (CFD) codes used to calculate prac- 
tical flows relevant to aerospace and aero-propulsion systems. 

2. Work Accomplished 

In the following, work on model development and model applications will be 
reported briefly. 

2.1 Model development 

Because of the limitations of computer speed and storage, DNS (direct numerical 
simulation) is limited to flows of moderate Reynolds numbers and simple geometries. 
For the present and the foreseeable future, turbulence modeling is the only viable 
approach for the calculation of complex turbulent flows of engineering interest. In 
turbulence modeling, the k-e eddy viscosity model is the most widely used model 
in engineering calculations. Many flows have been calculated by A: - e models. A 
set of model constants have emerged from these computational experiences, giving 
rise to what is commonly referred to as the standard k-e model^’^. The standard 
k — e model can be directly used for free shear flows. For wall bounded flows, it 
can be used in conjunction with wall functions which prescribe the flow fleld at the 
equilibrium log layer instead of at the wall. However universal wall functions do 
not exist in complex flows and it is thus necessary to develop a form of fc - e model 
equations which can be integrated down to the wall. Jones and Launder^ were the 
first to propose such a low Reynolds number k-e model for near wall turbulence, 
which was then followed by a number of similar k-e models. A critical evaluation 
of the pre-1985 models was made by Patel et al.^. More recently proposed models 
are found in Lang and Shih® . 

The existing k-e eddy viscosity models for near wall turbulence suffer from some 
well known deficiencies. (Some of the existing models may be free from one or two 
of these deficiencies.) First, a near wall pseudo-dissipation rate was introduced to 
remove the singularity in the dissipation rate equation at the wall. However, the 
definition of the near wall pseudo- dissipation rate was quite arbitrary. Second, the 
model constants were different from those of the standard k-e model, making 
the near wall models less capable of handhng flows containing both high Reynolds 
number turbulence and near wall turbulence, which is often the case for a real flow 
situation. Patel et al.^ hsted as the first criterion the abihty of the near wall models 



8 


Z. Yang 


to be able to predict turbulent free shear flows. Third, the variable ?/+ was used 
in the damping function of the eddy viscosity formula. Since the definition of 
involves Ut-, the friction velocity, any model containing can not be used in 
flows with separation and reattachment. In addition, may not be well defined 
for flows with complex geometries. 

All the k - e models for near turbulence reviewed in references 4 and 5 were 
constructed based on the standard k — e model, i.e., they took the standard k — e 
model form as the base form and added on it modifications due to the near wall 
effect. The standard k — e model form itself has the following deficiencies: 1) the 
model is not realizable in that physically unrealistic Reynolds stresses could be 
obtained from the model in certain flow situations; 2) contrary to the experimental 
findings, the predicted spreading rate for a round-jet is larger than that for a planar- 
jet, a phenomenon commonly known as the planar-jet/round-jet anomaly; 3) the 
inadequate predicted response to pressure gradient of turbulent boundary layer 
flows with adverse pressure gradients is also attributed to the model form of the 
standard k — e model, as pointed out by Wilcox®. 

The purpose of our research in the area of eddy viscosity modehng is to propose 
a. k — e eddy viscosity model which is free from the deficiencies mentioned above. In 
addition, we require the proposed model to be free from coordinate parameters in the 
model formulation. While coordinate information, distance to the wall for example, 
provides a convenient parameter to calibrate the model’s damping functions, its 
use makes the model not tensorially invariant and thus unable to be incorporated 
into computational fluid dynamics (CFD) codes with unstructured grids. The final 
product of the present research is aimed at a realizable, tensorially invariant, and 
Galilean invariant k — e eddy viscosity model which has the capabihty for both free 
shear flows and wall bounded flows. 

The components for such a model have been developed in the past few years. A 
variable formulation has been proposed, in which depends on the solutions of 
the mean flow field and the turbulence field. The formulation was derived based 
on the realizability constrains, rapid distortion theory (RDT), and the invariants 
theory. The Reynolds stresses from the k—e eddy viscosity model using the proposed 
formulation are always realizable in that no physically un-admissible Reynolds 
stresses will ever be produced. The variations of with flow situations are in the 
right direction as well, i.e., the value of is smaller than 0.09 in regions where it 
should be. Such a variation of with flow situations make the model calculation 
numerically more robust, particularly when the model is used to calculate flows with 
large strain rate variations, for example, the flow of shock-wave/turbulent boundary 
layer interactions. The detail of the proposed formulation and the corresponding 
model performance for benchmark free shear flows and wall bounded flows can be 
found in Shih et al.^ 

The dissipation rate equation in the standard k — e model has a singularity at the 
wall, where the turbulent kinetic energy is equal to zero. To remove this singularity, 
and yet without introducing an arbitrary pseudo-disspation rate, the k - e model 
was reformulated using velocity- scale and time-scale instead of the traditionally 
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used velocity-scaJe and length-scale. It was further argued that the turbulent time 
scale should be bounded from below by the Kolmogorov time scale which is always 
non-zero. With such a turbulence time scale, the resulting dissipation rate equation 
remains singularity-free everywhere in the flow field. The resulting dissipation rate 
equation is computationally very robust. The details of using Kolmogorov behaviors 
to aid the k-e eddy viscosity models for near wall turbulence can be found in Yang 
and Shih® and Shih and Lumley^. 

It is a challenging task to propose a fc - e eddy viscosity model for near wall 
turbulence which does not contain coordinate information. The difficulty lies in the 
lack of proper parameter (s) to calibrate the near wall effect. In their paper on the 
first k — e eddy viscosity model for near wall turbulence, Jones and Launder® used 
Rt, the turbulence Reynolds number, to calibrate the near wall effect. Rt is a field 
quantity, and does not depend on the coordinate system. Thus, a.k - e model with 
Rt is tensorially invariant. However, the performance of the resulting model is not 
very satisfactory, partly due to the fact that the near wall region for Rt is confined 
to y+ about 20 while the near region for the damping function extends to the 
lower end of the equilibrium log layer, which normally has a value of about 60. 
Other researchers have later used y, the distance to the wall, as a parameter to 
calibrate the near wall effect. Models using y in the damping functions are found to 
give a better performance due to the fact that y grows more gradually in the near 
wall region compared with Rt. However, any model with coordinate information 
will not be tensorially invariant. In addition, the definition of wall distance will be 
ambiguous in complex flows with complex geometries. 

Recently, Yang and Shih^° have introduced a new parameter R to calibrate the 
near wall effect. R is defined as R = ^ = and is related to the following 

two important physical parameters: the turbulent Reynolds number and the time 
scale ratio of the turbulence to the mean flow. R defined above is expressed in 
terms of the local field variables and is thus coordinate independent. In the near 
wall region, R is found to increase with in a gradual and monotonic manner. 
This property makes R an ideal candidate for constructing the damping function. 
In Yang and Shih^°, such a damping function was constructed to model the near 
wall effect. The proposed k-e model for near wall turbulence uses the standard 
k — e model as the corresponding high Reynolds number model and the resulting 
model was validated for turbulent channel flows and turbulent boundary layer flows 
with different pressure gradients. 

To adequately capture the effect of pressure gradient on turbulent boundary lay- 
ers, inhomogeneity effect was introduced in the dissipation rate equation. Since the 
dissipation rate represents the energy transfer from the large eddies to the small 
eddies in a turbulent flow, it is expected that the inhomogeneity in the large eddies 
would have an effect on the energy transfer rate. The previously used dissipation 
rate model is a homogeneous model in the sense that other than the diffusion term, 
the source terms remain the same for both homogeneous flows and inhomogeneous 
flows. In the present study, the effect of flow inhomogeneity is modeled directly 
as an extra production/destruction source term in the dissipation rate equation. 
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This extra source term is required to vanish when the flow is homogeneous. In the 
framework oik — e eddy viscosity model, and for the case of weak inhomogeneity, we 
characterize the inhomogeneity by Vfc and V5, which represent the inhomogeneity 
of the turbulent field and the inhomogeneity of the mean field, respectively. Invari- 
ants theory was used to deduce the possible forms of the source term. The resulting 
model constants in the source term were further calibrated against a turbulent 
boundary layer flow with favorable pressure gradients (the Herring and Norbury 
flow^^) and a turbulent boundary layer flow with adverse pressure gradients (the 
Bradshaw flow^^.) The resulting model was found to work well for other bound- 
ary layer flows with different pressure gradients. The model description and flow 
calculations can be found in Yang and Shih^^. 

The round-jet /planar-jet spreading rate anomaly is a long standing issue for tur- 
bulence modeling community. This anomaly is not only observed for two equation 
eddy viscosity models, it is also observed for second order closure models. This 
suggests that the dissipation rate equation needs to be modified in order to resolve 
this anomaly, as pointed out by other researchers. The present research confirms 
this conjecture and finds that either the extra source term in the dissipation rate 
equation due to the inhomogeneity effect reported in reference 13 or the new dis- 
sipation rate equation based on the dynamic equation for the fluctuating vorticity 
reported in reference 7 can resolve this round-jet /planar-jet spreading rate anomaly. 
The computational results are found in the corresponding references. 

Now, we have in place all the essential components to construct a realizable, ten- 
sorially invariant, and Galilean invariant A: — e eddy viscosity model which works well 
for both free shear flows and waU bounded flows. The variable formulation will 
ensure the realizability of the proposed model; the analysis of the Kolmogorov be- 
haviors of turbulence will make the model singularity-free without introducing any 
arbitrary pseudo-dissipation rate; a new parameter R has been found which possess 
the needed properties to calibrate the damping functions for near wall flows; the 
effects of pressure gradients on turbulent boundary layers are captured via an extra 
source term in the dissipation rate equation representing the contribution from flow 
inhomogeneity; the round-jet /planar-jet spreading rate anomaly can be resolved via 
the changes in the model dissipation rate equation, either by the contribution of the 
inhomogeneity effect or by the new dissipation rate equation based on the dynamic 
equation for the fluctuating vorticity. Now, what needs to be done is to put all the 
components together and test the resulting model against a large variety of flows 
of practical interest. This is a topic of our current research. 

2.2 Model applications 

The purpose of engineering turbulence modeling research is to use the proposed 
turbulence models to calculate flows of engineering interest. Thus, the ultimate test 
for our research on turbulence models hes in the abihty of the proposed models to 
calculate flows in aerospace and aero-propulsion systems. To this end, it is necessary 
to install the proposed models into a CFD code which serves as a numerical test-bed. 
In order to isolate the effect of turbulence model, the numerical test-bed should be 
a well established code with a well developed numerical scheme. The code should 
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be able to handle a wide range of flow situations and should be widely used in the 
aerospace and aero-propulsion community. 

In the present study, NPARC is chosen as the numerical test-bed for model ap- 
plications. NPARC is a general purpose computational fluid dynamics (CFD) code. 
It has a large user group in aerospace and aero-propulsion community. The code 
solves the Navier- Stokes equations in a conservation form in a general coordinate 
system. It can handle complex geometries and different types of boundary condi- 
tions. A description of the code and its numerical scheme is given by Cooper and 
Sirbaugh^^. Updates of code capabilities and a listing of bibhographies of using 
the NPARC code to calculate flows of practical interests can be found in the cur- 
rent version of the NPARC User Guide^^. The NPARC code is actively supported 
by NASA Lewis Research Center via the NPARC Alliance, a partnership formed 
between NASA Lewis Research Center and Air Force Arnold Engineering Devel- 
opment Center. The NPARC users are represented by the NPARC Association 
currently headed by Boeing Company and McDonnell Douglas Company. 

Two-equation k — e eddy viscosity turbulence models are incorporated into the 
NPARC code via a separate turbulence subprogram written by Dr. Zhu of CMOTT. 
(For a description of the philosophy and the structure of the subprogram, please see 
Zhu and Shih^®.) Different turbulence models have been tested against typical flows 
in propulsion systems, including, for example, the ejector nozzle flows, transonic 
diffuser flows, and the boat-tail nozzle flows. Details of flow calculations using the 
NPARC code and the turbulence subprogram are going to be presented by Yang et 
al.^^ in the upcoming 1995 AIAA Joint Propulsion Conference. 

By installing the advanced turbulence models into the NPARC code, NPARC’s 
turbulence modehng capability is enhanced. Turbulence model development for 
NPARC is CMOTT’s long term activity in supporting NASA’s High Speed Research 
(HSR) Program. Research in this area, together with model applications using the 
NPARC code, is described in more detail in this research brief in the article entitled 
“Turbulence Model Development for NPARC.” 

3. Future Plans 

1) Turbulence modeling: development and applications 

We will first finish the research on a realizable, tensorially invariant, and Galilean 
invariant k — e eddy viscosity model. The model is going to be validated for both 
simple shear flows, such as free shear flows and boundary layer flows, and the com- 
plex flows in aero-propulsion systems. Then, we will proceed with the development 
of Reynolds stress models for wall bounded turbulent flows. The final product in the 
latter area of research is a realizable, tensorially invariant, and Galilean invariant 
second order closure model. 

2) Transition modehng 

We will be working on the modehng and calculation of transitional flows over a 
low-pressure (LP) turbine blade. This research wih be part of NASA’s Low Pressure 
Turbine Flow Physics Program. For such flows, the transition is induced by the high 
turbulence level associated with the free-stream and the incoming wake, giving rise 
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to what is referred to as the bypass transition. The final product of this research 
will be a two equation model for transitional flow of bypass transition type. 
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Turbulence Modules for CFD Codes 

Jiang Zhu 


1. Motivation and Objectives 

It has been long recognized that there is a gap between turbulence model devel- 
opers and CFD users. The former mainly use simple flows to verify new modeling 
concepts and evaluate the resulting models, while the latter are usually reluctant to 
implement and test new, more advanced turbulence models unless and until they see 
such models showing good performance for a wide range of complex flow situations. 
Turbulence model equations also require special treatment to ensure numerical real- 
izability such as the positiveness of the turbulent kinetic energy and its dissipation 
rate. 

In order to bridge this gap, we develop turbulence modules for industry’s CFD 
codes. The modules are written in a self-contained manner so that the user can use 
any turbulence model in the modules without worrying about how it is implemented 
and solved. The input to a module is mean flow variables, boundary and geometric 
information which are to be provided by a mean flow solver, and the output of 
the module is the turbulent diffusivity and relevant turbulent source terms which 
are needed for the mean flow calculation. The interaction between the mean flow 
solver and the turbulence module will give the final turbulent flow solution. With 
the aid of the modules, we can also take the advantage of the well-estabhshed and 
sophisticated CFD codes to test turbulence models for a variety of complex flows 
which are intractable with the simple research codes. 

2. Work Accomplished 
2.1 The NPARC Code 

The NPARC code has been used extensively by the U.S. aerospace community 
to analyze propulsion flows. However until recently, only algebraic models such 
as the Thomas and Baldwin-Lomax have been available in NPARC for turbulent 
flow simulations. The Chien low Reynolds number k-e model with modifications 
for compressibility has been available in the 2D version but was not successfully 
installed in the 3D version. Another k-e model (the NPARC 1.0 k-c model) was 
installed in both the 2D and 3D versons but has not provided desirable accuracy 
and stability. 

In order to improve the NPARC’s capability to calculate turbulent flows, the 
two-equation turbulence model in the NAPRC code (both the 2D and 3D versions) 
was modified^ so that the model is based on the low Reynolds number k-e model of 
Chien and no longer on the NPARC 1.0 k-e model. Stability enhancements and a 
new inflow boundary condition for the turbulent quantities were also added to the 
k-e model. Comparisons of the NPARC solutions obtained using the previous and 
new models with experimental data indicated that the Chien k-e model installed in 
this work improves the capability of the NPARC to calculate propulsion flows. 
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2.2 A Turbulence Module for the FAST2D Code 

FAST2D is a computer program developed at the CFD group of Prof. Rodi^ 
for calculating two-dimensional, incompressible, elliptic flows with complex bound- 
aries. Here “two-dimensional” stands for plane or axisymmetric flows with or with- 
out swirhng. FAST2D is a finite-volume procedure and uses non-staggered grids 
and the cell-center arrangement, i.e., the flow variables are stored at the geometric 
centers of control volumes, with grid fines forming the faces of the control volumes. 
FAST2D is being used at CMOTT as a research code. The turbulence module, 
termed as FAST2DTS-1, is similar in structure to FAST2D. To reduce numeri- 
cal diffusion while maintaining necessary stability, the second-orde accurate HLPA 
scheme^ is used for the convective terms of the turbulent transport equations. This 
scheme is implemented such as to blend two schemes by a parameter A, with lim- 
iting values A=0 for the (first-order) upwind and A=1 for the second-order and 
bounded solution. Seven turbulence models have been built into the module: the 
standard k-e, CMOTT k-e, k-w, CMOTT Reynolds stress algebraic equation model. 
Launder- Sharma, Chien and Shih-Lumley, with the first four being high- and the 
last three being low-Reynolds number models. The module includes the four most 
commonly used boundary conditions - inflow plane, outflow plane, symmetry plane 
and solid wall. The solid wall boundary condition is formulated either by the stan- 
dard wall-function approach or by a low-Reynolds-number procedure, depending 
on the turbulence models used. The program also allows the user to introduce 
boundary conditions other than these four types. Care has been taken to make the 
module user-friendly. This is achieved by providing a few parameters defining the 
types and locations of boundaries, which allows the user to choose the boundary 
conditions given in the module by simply specifying these parameters. The system 
of algebraic difference equations is solved using the alternating direction TDMA of 
Thomas. 

2.3 A Turbulence Module for the NPARC 2D Code 

NPARC is a compressible flow code. Its 2D version is currently restricted to 
plane or axisymmetric flows without swirling. The turbulence module (NP2DTS- 
1)^ written for NPARC 2D has much in common with the above FAST2DTS-1. 
Only the following points are worthy of mention: (a) the cell-corner arrangement is 
used, i.e., the flow variables are stored at grid nodes rather than at the centers of 
control volumes; (b) contrary to NPARC and most other compressible codes, the 
non-delta form of equations is used which leads to simple linearization and is more 
effective to ensure the positiveness of the turbulent kinetic energy and its dissipation 
rate than the delta form; and (c) four turbulence models have been built into the 
module: Baldwin-Lomax, Chien, Shih-Lumley and CMOTT realizable models, the 
first being the algebraic eddy-viscosity model and the last three the low Reynolds 
number k-e models. The NPARC 2D code coupled with the module NP2DTS-1 have 
been applied® to a series of flows including the flow over a flat plate, in an ejector 
nozzle, in a transonic diffuser, and a boat-tail nozzle flow. Both NP2DTS-1 and 
FAST2DTS-1 were released at the 1994 Industry-Wide Workshop on Computational 
Turbulence Modeling, held at the Ohio Aerospace Institute. 


15 


Turbulence Modules for CFD Codes 

2.4 A Turbulence Module for the VSTAGE Code 

NASA LeRC currently has a very extensive research program on turbomachinery 
flow physics that is both experimental and computational. World class experimen- 
tal facilities exist to provide a multitude of experimental data for both component 
design needs and CFD comparisons. CFD turbomachinery tools are also well de- 
veloped after many years of testing and validation. The VSTAGE code, developed 
at Dr. Adamczyk’s group of the Lewis Research Academy, is one such tool. At 
present it is a proven CFD code which is being used by major U.S. engine com- 
panies in their design analysis. In a recent blind testing for a rotor compressor 
(designated Rotor 37) organized by the ASME/International Gas Turbine Confer- 
ence held at the Hague, 1994, the VSTAGE prediction turned out to be among the 
best of predictions submitted by a total of 12 groups. However, the best predictions 
for Rotor 37 were still less than satisfactory for the purpose of design, and this was 
mainly attributed to turbulence modehng. The VSTAGE is equipped only with 
the commonly used Baldwin-Lomax mixing length model. The limitations of such 
model are quite apparent from the predictions for Rotor 37 case. It is questionable 
whether improvements in predicting transonic rotor flows can be made using this 
model. Thus a turbulence module is developed for VSTAGE. The reasons of de- 
veloping this module are twofold: 1) to provide VSTAGE with more general and 
advanced turbulence models beyond the algebraic level; and 2) using VSTAGE as 
a common numerical platform to assess the performance of different existing and 
newly developed turbulence models for turbomachinery flows. The module has the 
same numerical framework as that of the NPARC module but is written exclusively 
for VSTAGE. As a result, it has the following unique features: (a) the cylindrical 
coordinate system (x, r, 6) is used as the absolute (flxed) reference frame, and cor- 
respondingly the cylindrical coordinate components of the flow velocity are used in 
the curvihnear transformation of equations; (b) the rotational effect is taken into ac- 
count by transforming the absolute reference frame to the relative (rotating) frame; 
(c) instead of the commonly used metrics such as the module uses vol- 

umes and areas of control volumes as the major geometric quantities, which is in 
lin e with the definition of variables in VSTAGE; and (d) the data transfer between 
VSATGE and the module is via the Fortran common blocks. Three high Reynolds 
number k-e models have been implemented, one standard and two CMOTT re- 
alizable models. The CMOTT models differ from each other in the dissipation 
rate equation, one using the standard equation and the other the newly developed 
equation®. The wall function approaches are used to formulate the sohd wall bound- 
ary conditions. The code validation is being under way, and the preliminary results 
for Rotor 37 test case look promising. 

2.5 Calculation of Wall-Bounded Complex Flows 

This study’’ concentrates on complex turbulent shear flows which are of great 
interest in propulsion systems. These flows are backward-facing step flows, con- 
fined coflowing jets, confined swirling coaxial jets, U-duct flows and diffuser flows. 
Most of these flows have complex structures. For example, the confined coflowing 
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jet combines several types of flow structures, such as the shear layer, jet, recircu- 
lation, separation and reattachment. Accurate prediction of these flows is of great 
importance for engine design in all its key elements. 

The turbulence model used in this study is a recently developed realizable Reynolds 
stress algebraic equation model®’® which is fundamentally different from the tradi- 
tional algebraic Reynolds stress models. The present model is developed using the 
invariance theory in continuum mechanics. This theory leads to a general constitu- 
tive relation for the Reynolds stress tensor in terms of the mean deformation rate 
tensor and the turbulent velocity and length scales characterized by the turbulent 
kinetic energy and its dissipation rate. Realizability is imposed directly on the 
constitutive relation for the Reynolds stresses to determine the coeflSicients in the 
relation. As a result, a realizable explicit expression for the Reynolds stresses is 
obtained for general three-dimensional turbulent flows. Some model constants were 
fine-tuned against a backward-facing step flow and then tested in other flows. 

The calculations were performed using the FAST2D code. Grid independent 
and low numerical diffusion solutions were obtained by using differencing schemes 
of second-order accuracy on suflSciently fine grids. The standard wall function 
approach was used for wall boundary conditions. Calculations using the standard 
k-€ (SKE) model were also carried out for the purpose of comparison. 

Diffuser Flows. Two conical diffuser flows were calculated, one with a 8° total 
ajigle (Trupp et ai.’s case) and the other 10° (Fraser’case). In both cases, the flows 
undergo strong adverse pressure gradients but remain attached. Although the flow 
conflguration looks simple, it is not easy to calculate this type of flow accurately, 
especially for the boundary layer quantities. Fig.l shows the variation of calculated 
and measured wall friction coefficient Cf with the axial distance x/Ro {Ro is the 
inlet duct radius). It is seen that the result of the present model is in good agreement 
with the experimental data, while the SKE model overpredicts Cf along almost the 
entire length of the diffuser. The calculated and measured displacement thickness 
6* are compared in Fig.2. The comparison shows that the SKE model gives a good 
prediction in the upstream region, but deviates significantly from the experiment 
downstream; the present model prediction is good in the whole region. Fig.3 shows 
the comparison of calculated and measured shape factor H. This is the case in 
which the worst agreement with the measurement has been found for both models. 
Nevertheless, the present model still performs considerably better than does the 
SKE model. 

U-Duct Flow. This case is the experiment of Monson et al. (1990) conducted 
in a 180° planar turnaround duct. It features flow with large streamline curvature. 
Calculations were compared to the experiment taken at a flow Reynolds number 
of 10®. Fig.4 shows the streamlines computed with the present model. A small 
separation region is found at the bend exit. However, the SKE model did not 
predict the flow separation. Fig.5 shows the comparison of calculated and measured 
Cf along the inner wall. The bend is located between 21.7<s/H<24.8. Both models 
are seen to behave in the same manner and produce large discrepancies in the bend 
region. The reason for this may partially due to the use of the wall function which 
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does not respond to the severe pressure gradient. 

Backward- Facing Step Flows. Two backward facing step flows, measured by 
Driver and Seegmiller (1985) and Kim et al. (1978), were calculated. The former 
(DS case) has a smaller and the latter (KKJ case) a larger step expansion. The com- 
puted and measured reattachment points are compared in Table 1. The calculated 
reattachment point from the present model agrees well with the experiments. Fig.6 
shows the comparison of the computed and the measured static pressure coefficient 
Cp along the bottom wall. The SKE model is seen to predict a premature pressure 
rise, which is consistent with its underprediction of the reattachment length, while 
the present model captures the pressure rise quite well. Fig. 7 shows the compar- 
isons of predicted and measured turbulent stresses uu, vv and uv at the location 
x=2 which is in the recirculation region. In the KKJ-case, no reliable experimental 
data exist for the turbulent stresses due to the unsteadiness of the flow. However, 
the experimental data of the DS-case is considered more reliable because of the 
smaller unsteadiness of the flow. As compared with the results of the SKE model 
in Fig.7, it is seen that the anisotropic terms in the present model increase uu and 
decrease W, leading to significant improvements in both uu and vv except in the 
near-wall region. On the other hand, the anisotropic terms have little impact on 
uv. The improvement obtained by the present model for uv is mainly due to the 
reduction in Cp by strain rate. 


Table 1. Comparison of the reattachment points 


Case 

measurement 

SKE 

PRESENT 

DS 

6.26 

4.99 

5.82 

KKJ 

7± 0.5 

6.35 

7.35 


Confined Jets. The general features of confined jets are sketched in Fig.8. At 
the entrance, two uniform flows, a jet of larger velocity and an ambient stream of 
smaller velocity, are discharged into a cylindrical duct of diameter Do- The inlet flow 
conditions can be characterized by the Craya-Curtet number Ct- The experiment 
shows that recirculation occurs when Ct <0.96. For a given geometry, recircula- 
tion as weU as adverse pressure gradients can be intensified by reducing the value 
of Ct at the entrance. The separation and reattachment points of the predicted 
recirculation bubbles are compared with the experimental data in Fig.9. The ex- 
periment indicated that as Ct decreases, the separation point moves upstream while 
the reattachment point remains practically unchanged. The present model captures 
this feature well and predicts both the separation and reattachment points much 
better than does the SKE model. The variation of the pressure coefficient Cp along 
the duct wall is shown in Fig. 10. The pressure distribution is governed by the jet 
entrainment as well as the contraction and expansion of the flow caused by the 
recirculation bubble. The decrease in the ambient velocity induced by the entrain- 
ment gives rise to an adverse pressure gradient, while the contraction of streamlines 
produces the opposite effect. These two mechanisms interact more intensely with 
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each other as Ct decreases and cause the pressure to vary httle in the region up- 
stream of the center of the recirculation bubble. However, in the downstream part 
of the recirculation bubble, the deceleration of the flow sets up an adverse pressure 
gradient, the slope of which becomes steeper as Ct decreases. Therefore, the ability 
to capture the location of the recirculation center will have a direct impact on the 
prediction of the pressure. Regarding the comparison between predictions and ex- 
periments, it is seen that although both models predict practically the same total 
pressure rises which are in excellent agreement with the measurements, the present 
model captures the steep pressure gradients better than does the SKE model for all 
of the Ct values. 

Confined Swirling Coaxial Jets. This is the case experimentally studied 
by Roback and Johnson (1983). Fig.ll shows the general features of the flow. 
At the inlet, an inner jet and an annular jet are ejected into an enlarged duct. 
Besides an annular recirculation bubble due to sudden expansion of the duct, a 
centerline recirculation bubble is created by flow swirling. Fig. 12 compares the 
calculation of the centerline velocity with the experiment. The negative velocity 
indicates the central recirculation. It is seen that both models predict the strength 
of central recirculation and the front stagnation point quite well, but the present 
model predicts the rear stagnation point much better than does the SKE model. 
Fig. 13 shows the comparison of calculated and measured mean velocity profiles at 
x=5.1cm. Both models give reasonably good profiles which are within experimental 
scatter, except for the peak values of the axial and radial velocities. Both models 
have been found to give nearly the same results in the downstream region, which 
can also be seen from Fig.l2. 

These comparisons show that the present realizable Reynolds stress algebraic 
equation model significantly improves the predictive capability of k-e equation based 
models, expecially for flows involving massive separations or strong shear layers. In 
these situations, the standard eddy viscosity model overpredicts the eddy viscosity 
and, hence, fails to accurately predict wall shear stress, separation, recirculation, 
etc. We find that the success of the present model in modeling the above men- 
tioned complex flows is largely due to its eff^ective eddy viscosity formulation which 
accounts for the eflFect of mean shear rates. According to the present model, the 
effective eddy viscosity will be significantly reduced by the mean strain rate and 
maintained at a correct level to mimic the complex flow structures. 

2.6 A New k-e Eddy Viscosity Model 

The model® has a new dissipation rate equation and a new realizable eddy vis- 
cosity formulation. The former is based on the dynamic equation for fluctuating 
vorticity and the latter is derived from the realizability analysis. The model contains 
the effect of mean rotation on turbulent stresses, and its dissipation rate equation 
has an always positive production term, which is expected to enhance the numerical 
stability in turbulent flow calculations, especially when using more advanced closure 
schemes such as second-order closures. 

Comparisons with experimental data show that the model performs better than 
the standard k-e model for almost all the flows tested®. The performance of the 
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model for complicated recirculating flows is demonstrated through calculations for 
two backward-facing step flows, one (DS-case) with smaller and the other (KKJ- 
case) with larger step height, both of which have been extensively used in the 
literature to benchmark calculations of separated flows. Sufficiently fine grids, with 
201x109 points in the DS-case and 199x91 points in the KKJ-case, were used to 
establish numerical credibility of the solutions. The computational domain had a 
length of 50 step heights, one fifth of which was placed upstream of the step. The ex- 
perimental data were used to specify the inflow conditions, the fully-developed flow 
conditions were imposed at the outflow boundary, and the standard wall function 
approach was used to bridge the viscous sublayer near the wall. The compari- 
son with the standard k-e model and the experiments is shown in Table 2 for the 
reattachment length. 


Table 2. Comparison of the reattachment point locations 


Case 

measurement 

standard model 

present model 

DS 

6.26 

4.99 

6.02 

KKJ 

7± 0.5 

6.35 

7.5 


The comparison of the size of the separation bubble, the skin friction pressure 
coefllcients shows that the overall performance of the present model is better than 
that of the standard model. 

2.7 A New Reynolds Stress Algebraic Equation Model 

A new Reynolds stress algebraic equation model has been developed using a 
truncated constitutive relation^°. This model differs from the previous one® mainly 
in: (a) it has a simpler quadratic form; (b) it shows the proper lack of a rotation 
effect on the isotropic turbulence, satisfying the rapid distortion theory; and (c) it 
is fully realizable, i.e., it ensures both the positivity of the normal Reynolds stresses 
and the Schwarz’ inequality between turbulent velocity correlations. The model 
were first calibrated using well-studied basic flows such as homogenous shear flow 
and the surface flow in the inertial sublayer and then applied to complex flows 
including the separated flow over a backward-facing step and the flow in a confined 
jet. 

Basic flows. Two basic cases were considered; they are Tavoularis and Corrsin’s 
(1981) homogeneous shear flow and the direct numerical simulation of channel flow 
(Kim, 1990). The model gives 612 = —0.156, bu = —^22 = 0.123 for the homo- 
geneous shear flow at Ui^ikje = 6.08 and gives b \2 = —0.122, 6 n = —622 = 0-14 
for the channel flow in the inertial sublayer at Ui^^kje = 3.3. These results show 
that the present model gives reasonable anisotropy of Reynolds stresses for both 
the homogeneous shear flow and the boundary layer flow compared to the standard 
k-€ eddy viscosity model which gives 5n = 622 = 0 for both the flows and gives 
bi 2 = —0.273 for the homogeneous shear flow and 612 = —0.149 for the boundary 
layer flow. Detailed comparisons with the experimental and DNS data are shown in 
Table 3 for the homogeneous shear flow and in Table 4 and Fig. 14 for the channel 
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flow. 


Table 3. Anisotropy in the homogeneous shear flow 



experiment 

standard 

present 

6i2 

-0.142 

-0.273 

-0.156 

hi 

0.202 

0 . 

0.123 

622 

-0.145 

0 . 

-0.123 


Table 4. Anisotropy in the channel flow 



DNS data 

standard 

present 

612 

-0.145 

-0.149 

- 0.122 

hi 

0.175 

0 . 

0.14 

622 

-0.145 

0 . 

-0.14 


Backward-facing step flows. Figs.l5(a) and (b) show the friction coefficient 
C/ at the bottom wall calculated with the SKE model and the present model; also 
included in fig. 15(a) are the experimental data for the DS-case, but no such data 
are available for the KKJ-case. It can be seen that the grid refinement does produce 
some differences for the results of the present model, more noticeable in the KKJ- 
case, and this is also the case for the SKE results. This indicates that the solutions 
obtained on the coarse grids are not sufficiently close to the grid-independent stage. 
Recently, Thangam and Hur (1991) have conducted a highly-resolved calculation 
for the KKJ-case. They have found that quadrupling a 166x73 grid leads to only 
a minimal improvement. Therefore, the present results on the fine grids can be 
considered as grid-independent. For the DS-case, the fine grid computations with 
the SKE model and present model required 703 and 691 iterations, and took ap- 
proximately 7.1 and 9 minutes of CPU time on the Cray YMP computer. In the 
following, only the fine grid results are presented. 

The wall friction coefficient Cf is a parameter that is very sensitive to the near- 
wall turbulence modeling. It is C/ that the various low Reynolds number k-e models 
tested predict much worse than those using wall functions. However, the influence 
of the near-wall turbulence modeling is mainly restricted to the near-wall regions. 
It is seen from fig. 15(a) that both the SKE model and the present model largely 
underpredict the negative peak of Cf, pointing to limited accuracy of the wall 
function approach in the recirculation region. 

The computed and measured reattachment points are compared in Table 5. They 
are determined in the calculation from the point where Cf goes to zero. The 
reattachment point is a critical parameter which has often been used to assess the 
overall performance of turbulence models as well as numerical procedures. Table 5 
clearly demonstrates the significant improvement obtained with the present model. 
It is important to mention that this improvement is mainly due to the behavior of 
Cfj, in the present model, and that the anisotropic behavior of the turbulent stresses 
only makes a marginal contribution to it. 
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Table 5. Comparison of the reattachment points 


Case 

measurement 

SKE 

PRESENT 

DS 

6.26 

4.99 

5.80 

KKJ 

7± 0.5 

6.35 

7.27 


Figs. 16(a) and (b) show the comparison of computed and measured static pressure 
coefficients Cp along the bottom wall. In both cases, the SKE model is seen to 
predict premature pressure rises which is consistent with its underprediction of the 
reattachment lengths. 

Finally, the comparisons of predicted and measured turbulent stresses uu, vv 
and uv are shown in figs.17 and 18 at various x-locations. In the KKJ-case, no 
experimental data for the turbulent stresses are available in the recirculation region, 
and the reattachment point was found in the experiment to move forward and 
backward continuously around seven step heights downstream of the step, leaving 
an uncertainty of ±0.5 step height for the reattachment length. This also points to 
some uncertainty in the measured turbulent quantities in the recovery region. On 
the other hand, the experimental data in the DS-case should be considered more 
reliable because of the smaller uncertainty of the reattachment location, indicating 
a smaller unsteadiness of the flow. The SKE model gives unrealistic results about 
normal Reynolds stresses: uu > titi at all the locations. In contrast, the present 
model gives at least qualitatively correct results due to the non-linear terms which 
increase uu while decreasing uu, leading to an overall improvement in both uu and 
vv results. 

Confined Jets. The predicted axial mean velocity profiles at two C* numbers 
are shown and compared with the experimental data in fig. 19, where R and Um 
are the radius of the cylinder and the sectional mean velocity, respectively. Both 
models are seen to predict very well the upstream evolution of the flow. As for the 
downstream development, the results of the present model remain in good agreement 
with experiments while the SKE model underpredicts the centerline velocity decay 
at all Ct numbers. The variation of the pressure coefficient Cp along the duct wall 
is shown in fig.20. The pressure distribution is governed by the jet entrainment 
as well as the contraction and expansion of the flow caused by the recirculation 
bubble. The decrease in the ambient velocity induced by the entrainment gives rise 
to an adverse pressure gradient, while the contraction of streamlines produces the 
opposite effect. These two mechanisms interact more intensely with each other as Ct 
decreases and cause the pressure to vary little in the region upstream of the center 
of the recirculation bubble. However, in the downstream part of the recirculation 
bubble, the deceleration of the flow sets up an adverse pressure gradient, the slope 
of which becomes steeper as C< decreases. Therefore, the ability to capture the 
location of the recirculation center will have a direct impact on the prediction of 
the pressure. Regarding the comparison between predictions and experiments, it is 
seen that although both models predict the same total pressure rises which are in 
excellent agreement with the measurements, the present model captures the pressure 
distribution much better than does the SKE model for all the Ct values. 




22 


J. Zhu 


3. Future Plans 

1) To develop a 3D version of the turbulence module for the NPARC code; 

2) To develop and validate the turbulence module for the VSTAGE code; 

3) To test turbulence models developed at the CMOTT and elsewhere. 
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Fig.l4 Direct comparison with DNS data of channel flow 
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Fig.l5 Friction coefficient along the bottom wall. 

: present model, fine grid; : present model, coarse grid; 

: SKE, fine grid; •: experiment 





Fig. 16 Pressure coefficient along the bottom wall. 
: present model; : SKE; •: experiment 
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Fig. 17 Turbulent stress profiles in DS-case. 
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Fig. 19 Axial mean velocity profiles. 
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Fig.20 Pressure coefficient along the duct wall. 
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Eddy Diffusivity Model for 
Turbulent Heat Transfer 

Aamir Shabbir 

1. Motivation and Objectives 

Commonly, the turbulent heat transfer is calculated by assuming a constant value 
of the turbulent Prandtl number. A mixing length turbulence model or a two 
equation k - e model is usually used in such calculations to calculate the turbulent 
eddy viscosity. Obviously this approach has hmitations in situations where the 
turbulent Prandtl number is not constant. Realizing this, some recent studies have 
proposed new eddy diffusivity models in which the transport equations for the scalar 
variance and its dissipation rate are solved to calculate the thermal diffusivity. 
These include, for instance the work of Nagano and Kim^ and Youssef, Nagano 
and Tagawa^ Hattouri, Nagano, and Tagawa^. The last study presents an updated 
and refined version of Nagano and Kim model. There are several other studies 
reported in the literature where the near wall modehng issues are discussed. In 
this report only the high Reynolds number version of the models is discussed, and 
therefore, reader should consult other papers which deal with this issue (for instance 
the review paper by So^). Recently Schwab and Lakshminarayana^ also propsed 
a two equation model for scalar field in which a scalar time scale equation is used 
instead of the scalar dissipation rate. They sucessfuUy apphed their model to a host 
of homogeneous flows. 

In last year’s research briefs I descirbed my project aimed at developing a new 
two equation eddy diflFusivity model for calculating turbulent heat transfer. In this 
model transport equations for the scalar variance and its dissipation rate are used to 
calculate the turbulent eddy diffusivity. As was pointed out then, the present effort 
differs from the other work in the following two respects. (1) In the above cited 
works, the extension of the scalar dissipation rate equation is based upon the work 
of Newman et al^ who developed its production/destruction mechanisms analogous 
to those of the mechanical dissipation rate equation. The model equation proposed 
in the present study is based on the exact transport equation for scalar dissipation 
and, its production/destruction mechanisms differ from those proposed in the other 
studies. (2) The model coefficient in the the scalar flux constitutive relation used 
in the present study is not a constant but is a function of the local invariants. 

2. Work Accomplished 

2.1 A New eddy viscostiy fc — e model for High Reynolds Number Flows 

Some of the effort was directed in developing and assessing the performance of 
a new two equation k-epsilon model in a joint effort with the other CMOTT re- 
searchers. The details of the model development and application can be found in 
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the paper by Shih et al J The model performs reasonably for a host of benchmark 
flows which includes: rotating homogeneous shear flow; various flat plate boundary 
layers; mixing layers; and jets. We, in this section, only show the apphcation of the 
model to the rotating homogeneous sheax flow, for which a large eddy simulation 
was carried out by Bardina et al}. We will compare the results from the present k-e 
model as well as the standard k-e model with the LES data for four different cases 
of f2/S (which are O/5=0.0, f2/5’=-0.50, 11/5=0.25, and 11/5=0.50). The initial 
conditions in all these cases correspond to isotropic turbulence and eoZ-S'fco = 0.296. 
Figure 1 (a) compares the evolution of turbulence kinetic energy, normalized by its 
initial value with the non-dimensional time St for the case of H/5 = 0.0. For 
this case both the present and the standard k — e (denoted by ske hereafter) models 
show the trends exhibited by LES, with the present model closer to the LES data. 
Figure 1 (b) shows the comparisons for the case 11/5 = 0.25. The LES shows that 
the growth rate of the turbulence kinetic energy is increased over the no rotation 
rate case. The present model is able to pick up this trend while the ske model does 
not. Figures 1 (c) and 1 (d) compare the evolution of turbulence kinetic energy for 
two more cases of H/5 = 0.5 and 11/5 = —0.5. For the first of these cases the LES 
shows that the growth rate of the turbulence kinetic energy decreases over the no 
rotation rate case. The present model is able to pick up this trend and although 
the agreement between the present model and the LES is not as good as it is for the 
other cases, it still is a lot better than the ske model. For the case of H/5 = —0.5 
the ske model does not show the effect of rotation on turbulence as it gives the same 
growth rate of turbulence kinetic energy as it did for the no rotation case, a result 
which is already known. On the other hand the present model is in reasonable 
agreement with the LES data as it shows the decay of the turbulence kinetic energy 
with time. 

2.2 A New eddy diffusivity 9^—ee model for High Reynolds Number Flows 

The details of the development of the scalar (or thermal) dissipation rate equation 
and the scalar flux can be found in the last year’s research briefs. Since then the 
model has been extended to inhomogeneous flows by modeling the diffusion term 
in the disspiation rate equation. A simple gradient diffusion type model is used for 
the diffusion term {uje$),j. With this the scalar dissipation rate equation is given 
by 


),j +CeieeS + ^ — , (1) 

<^<i> y/Pr [k -t- y/ue) 

where the model constant is to be determined. To do so we consider the log-law 
region of a wall, where the turbulent scalar (or heat) flux is constant. In this region 
the advection term can be ignored and equation (l) reduces to 


d fOLT dee^ ^ ^ c , ^ 

CeiegS -I- Cg2 - 7^ $ — Cqzti~. — 


dy dy 


y/Pr 


{k + y/ve) 


=0 


( 2 ) 


Eddy Diffusivity Model for Turbulent Heat Transfer 


35 


Using the log-law relations, and the fact that in the log-layer production equals 
the dissipation rate, the above equation gives the following relation for the model 
const 3ft ^ (p 


- Cn - 

where the Von Karman constant K is 0.41 and Ke is 0.482^. We know that in 
the log-layer = 0.09 and Cx = 0.1. With these values we obtain = 1.8 
from (3). The expression for the model constant C\ used previoulsy was somewhat 
cumbersome. We have now considerably simplified it and is given by equation (7). 
The model for the scalar (or thermal) field can now be summarized by the following 
equations. 
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Sr = x/Q^i^ S = ^2S;j^j, v = Skje, ^ f 
Cei = Cl - 0.33, Cs 2 = 0.63, C^ = C 2 - 14- r, at = 1.0, 

2.3 Application of Model 

Since the expression for C\ has been simplified than what was reported previously, 
it is instructive to re-calculate the fiows which were calculated previously. Figures 
1-7 show the evolution of the temperature variance and its dissipation rate for the 
seven different cases of Sirivat and Warhaft^° experiment. We note that the present 
model is in reasonable agreement with the experiment, and performs better than 
the Hattouri, Nagano and Tagawa model (HNT). The results on the overall are 
same as what was reported previously. Same is true for the homogeneous shear flow 
with a constant temperature gradient, for which the results are shown in figure 8. 

The model was also applied to the flat plate boundary layer whose surface is 
heated to a constant temperature. The experiment of Gibson et al.^ provides one 
test case for such a flow. This flow constitutes an initial value problem in the stream- 
wise direction and a boundary value in the cross-stream direction. It was calcualted 
using a quasi-implicit finite difference scheme. Since only the high Reynolds number 
form of the present model is developed so far, it is used here in conjunction with the 
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wall functions. For the k and e quantities the standard wall functions were used. 
The wall functions used for the thermal field are summarized in Appendix B. (It 
is found that with these wall functions the log-layer temperature profile agrees well 
with the data but the Stanton number is underestimated by about 7%. Work is 
underway to address this aspect of the wall functions.) The equations set consisted 
of: mean momentum and mean energy equations; turbulence kinetic energy and its 
dissipation equation; temperature variance and its dissipation equation. 

The results for the mean temperature are shown in figure 9. We note that the 
agreement with the experiment is reasonable. For comparisons purposes results 
are also shown for the HNT model which agrees with the experiment very well. 
Note that the HNT model was not used with the wall functions but was instead 
integrated to the wall. 

3. Future Plans 

The present model will be extended for integration to the wall. This will neces- 
sitate the formulation of damping functions for the model. However, due to the 
changing research environment at NASA, most of the next year will be spent on 
application of two- equation models to the turbomachinery flows. 
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Appendix A 

The k — e model used in the present study is that of Shih et alJ and is summarized 
below. 
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Note that the Tlij is the mean rotation rate as viewed in the rotating reference 
frame and uJk is the angular velocity at which the reference frame is rotating. The 
constants in the e equation are: Ci = 0.42 and C 2 = 1.9. 
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Appendix B 

This appendix gives the wall functions which were used to obtain the boundary 
conditions for the thermal field (i.e the values of 0 , and ee at the first grid point 
off the wall). Note that the standard wall functions were used for the velocity field 
and, therefore, are not listed here. 

For the mean temperature the standard log-law for a flat plate with constant 
surface temperature was used as a boundary condition (see e.g. Launder^^) and is 
given cis 


where A© is the difference between the wall temperature and the local temperature. 
From the experiment of Gibson et al.^ Ke = 0.482 and Ce = 3.8. For the turbulent 
part of the temperature field the local equilibrium assumption was invoked (i.e. 
Pe = € 6 ») to obtain the following relation for the boundary conditions on 

-=Cy^jC^ (B2) 


Once 02 is known the boundary condition on ee can be readily obtained by using 
the definition of the time scale ratio. 




(B3) 


Note that for the log-layer the time scale ratio r is 2.0; = 0.09; and C\ = 0.1. 
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Figure 1. Evolution of the normalized temperature variance and thermal dissipation 
rate in the experiment of Sirivat and Warhaft (1983). The temperature variance 
is normalized as = 6^J6^q and the thermal dissipation is normalized as 6^ 

— . For this case Wo = 0.0128 uq = 0.145515m/s, and lo = 0.011937m. 



Figure 2. Evolution of the normalized temperature variance and thermal dissipation 
rate in the experiment of Sirivat and Warhaft (1983). The temperature variance 
is normalized as = W/Wq and the thermal dissipation is normalized as = 
. For this caseFo = 0.002287 uo = 0.145515m/s, and k = 0.011937m. 
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Figure 3. Evolution of the normalized temperature variance and thermal dissipation 
rate in the experiment of Sirivat and Warhaft (1983). The temperature variance 
is normalized as 6?^ = and the thermal dissipation is normalized as = 

— . For this caseFo = 0.001705 uq = 0.145515m/s, and Iq = 0.011937m. 



X/M 


Figure 4. Evolution of the normalized temperature variance and thermal dissipation 
rate in the experiment of Sirivat and Warhaft (1983). The temperature variance 
is normalized as = 9^ f 6^ q and the thermal dissipation is normalized as = 
— . For this case = 0.009059 uq = 0.074701m/s, and Iq = 0.011071m. 
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Figure 5 Evolution of the normalized temperature variance and thermal dissipation 
rate in the experiment of ^vat and Warhaft (1983). The temperature variance 

is normalized as ^ = ¥/e'^o and the thermal dissipation is normalized as e’g = 
. For this case = 0.000924 tio = 0.074701m/s, and lo = 0.011071m. 



Fiffure 6 Evolution of the normalized temperature variance and thermal dissi- 
pation rate in the experiment of_Sinvat and Warhaft (1983). The temperature 

variance is normalized as 9*' = and the thermal dissipation is normalized 

For this case Bh = 0.0004471 “C=, no = 0.074701m/s, and 


as e 


■I _ 


Iq = 0.011071m. 
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Figure 7. Evolution of the normalized temperature variance and thermal dissi- 
pation rate in the experiment of Sirivat and Warhaft (1983). The temperature 
variance is normalized as and the thermal dissipation is normalized 

as ea = — — . For this case = 0.0004955 uq = 0.074701m/s, and 

io = 0.011071m. 



Figure 8a. Evolution of the normalized temperature variance in the experiment 
of Tavoularis and Corrsin (1981). The temperature variance is normalized with 
^^«/ = 0.01 
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Figure 8b. Evolution of the thermal dissipation rate in the experiment of Tavoularis 
and Corrsin (1981). 
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Figure 9. Mean temperature profile in wall units for the constant temperature flat 
plate boundary layer for the experiment of Gibson et al? . 
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Computations of Compressible 
Turbulent Benchmark Flows 

William W. Liou 

1. Motivation and Objective 

The development of advanced propulsion systems for high speed aerospace vehi- 
cles will require accurate computational models of turbulence that can be used in 
the CFD calculation of individual flow components. Before a new model is used in 
the calculation of complex flows in propulsion systems, the model has to be assessed 
in simple benchmark flows that not only contain the essential physical mechanisms 
at work in engine component flows, but are also well-documented. An immediate 
objective of the research activity described here is to identify these compressible 
benchmark flows and performed a prehminary assessment of models developed lo- 
cally at CMOTT. The long term objective is to develop second-order closure models 
for compressible flows. 

2. Work Accomplished 

2.1 Model Validation in Compressible Flows 

Six diflterent types of two-dimensional flow fields typically encountered in aero- 
propulsion systems were selected. They are: (1) compressible free shear layers; (2) 
compressible boundary layers; (3) transonic bump flows; (4) supersonic ramp flows; 
(5) oblique shock- wave/turbulent boundary-layer interactions; (6) can combustor. 
Sketches of the flows are shown in Fig. 1. 

(1) The compressible free shear layer is an integral element in engine component 
flows, such as combustor and nozzle flows. In many instances, the turbulence mixing 
in the shear layers is a major factor in the operating efficiency of an engine compo- 
nent. Experiments seem to support a state of self-preserving in a fully developed 
shear flow, which is describable by local characteristic scales. It is also observed 
in experiments that the spreading rate of a compressible turbulent shear layer is 
far smaller than an incompressible shear layer with the same velocity and density 
ratios. The reduced turbulent mixing influences significantly the combustion effi- 
ciency and the noise emitted from high speed jet flows. Therefore, it is important 
that turbulence models can predict correctly the compressible free shear flow; in 
particular, its spreading rate. 

(2) To evaluate model performance in wall bounded compressible flows, models 
will be assessed in fully developed compressible turbulent boundary layers. It is gen- 
erally accepted that for free stream Mach numbers less than about 5, the turbulent 
structure in compressible boundary layers is nearly the same as in the corresponding 
incompressible (constant density) flows. This is Morkovin’s hypothesis. Therefore, 
with properly scaling, models that gives good results in incompressible flow may 
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also perform well in compressible boundary layer with free steam Mach numbers 
less than 5. 

The above two cases provide tests for models in a “clean” environment. The 
following four cases involve flows with either shock- waves or chemical reactions. 

(3) Case 3 resembles the transonic flow over an airfoil, which affects the per- 
formance of engine components such as compressors and turbines. The strong 
viscous/inviscid interactions in the flow, through a change in the effective body- 
shape, pose a critical test for turbulence closure models. Two flows were chosen. 
The first one corresponds to the experimental set-up of a flow studied by Delery 
et. al The boundary layer developed on the bump mounted on a wind-tunnel 
wall remains attached throughout the interaction. The flow speed becomes super- 
sonic atop the bump. The wind-tunnel is choked. This flow was also chosen in a 
recent EUROVAL^ eflFort to evaluate the response of a turbulent boundary layer 
developing on an airfoil to shock-waves. There, the case was designated as ON- 
ERA BUMP A. The second configuration corresponds to an experiment by Bachalo 
and Johnson^, where a toroidal bump is mounted on a cylinder. The axisymmetric 
arrangement eliminates the three-dimensional effects due to the boundary layer on 
the side walls of the wind tunnel. The boundary layer separates downstream of the 
shock wave. The axisymmetric configuration has made this flow particularly suit- 
able for the study of a flow with separation resulting from strong shock/turbulent 
boundary-layer interactions. 

(4) Cases 4 and 5 involve interactions between shock-wave and boundary layer 
in supersonic flows. The compression ramp flow is a classical problem that en- 
compasses many complex phenomena and is of great practical importance in the 
design of engine components, such as inlets and turbomachinery. The four ramp 
flows selected were investigated in a series of work by Settles and co- workers^. The 
free stream Mach numbers are less than 3 and, depending upon the ramp angle, 
flow separation may occurs at the corner of the ramp. The flows were proposed for 
model testing in shock wave/turbulent boundary-layer interactions in the 1980-1981 
AFSOR-Stanford Conference on complex turbulent flows. 

(5) Case 5 includes another fundamentally important flow that is critical to the 
performance of inlets and turbomachinery: the impingement of an oblique incident 
shock-wave on a wall®. A reflected shock is formed to turn the downstream flow 
parallel to the wall. Again, if the shocks were sufficiently strong, the turbulent 
boundary layer may separate. The challenge for a turbulence model is to predict 
the incipient separation of the boundary layer and the subsequent separated flow. 

(6) The flow in a engine combustor typically involves highly turbulent reacting 
regions. Therefore, a robust turbulence model is just as important as a quality 
chemistry model in the design of a high efficiency and low pollutant emission com- 
bustor. An axisymmetric can combustor configuration® tested in UC Irvine was 
examined. The experimental set-up includes an air-swirler and gaseous fuel injec- 
tion of propane. The fuel and air are not premixed before entering the combustor. 
This configuration is representative of a gas turbine engine combustor. 

These test cases are well-documented and should provide good data base for 
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objective evaluations of turbulence models. For instance, some experiments have 
reported results obtained using different measurement techniques, i.e. hot wire vs. 
LDV. In fact, some flows were also selected by previous large-scale model evaluation 
efforts. 

It is more efficiently to calculate the free shear layer and the flat plate boundary 
layer with a boundary-layer equation solver. When it is necessary to use a Navier- 
Stokes solver, the model solution of a flat plate boundary layer is compared with the 
solution obtained from the boundary-layer equation solver to validate the model im- 
plementation in the Navier-Stokes solver. The COMTUR code‘s is currently used in 
the solution of compressible Favre-average Navier-Stokes equations. The FAST2D 
code® is used for solving the variable density form of the Navier-Stokes equations. 
The chemical reactions were modeled by using a CMOTT scalar probability density 
function(PDF) code, LPDF2D. 

In the following, a few sample results of this model assessment effort are shown. 
More details can be found in NASA TM reports that are in preparation. 

Fig. 2 shows the skin friction and the wall pressure distributions for a supersonic 
flow over a 24 ramp^. The incoming free stream Mach number is 2.8. and the unit 
Reynolds number is 6.3xl0'’’/m. No compressibility corrections were used. The 
results show that the standard k — e model and the multiple-scale model^ give 
reasonable predictions for the flow recovery and the flow separation. RNG modeP® 
and the model predict a larger separation region than the measurement. 

Note that only the hnear part of the constitutive relation of the S&Z model was 
used in this calculation. 

Fig. 3 shows the variations of the skin friction coefficient and the wall pressure for 
an oblique shock/turbulent boundary-layer interaction. The incoming free stream 
Mach number is 2.89 and the unit Reynolds number is 5.73 x The impinging 

shock is strong enough to cause the boundary layer to separate. KEl modeP^-^®, 
which is a modified low Reynolds number model of Shih and Lumley^^, performs 
better than the Chien’s modeP®. The high Reynolds number KE2^® model and the 
standard k - e model gives reasonable predictions for both quantities. 

The predicted and measured wall pressure distributions for the ONERA BUMP 
A are shown in Fig. 4. KE2^® model give a slightly better predication downstream 
of the shock wave than the other models. All the model predictions asymptote to a 
value higher than the measurement in the fully developed region. This may be due 
to three-dimensional effects, such as the boundary layer on the wind-tunnel walls. 

The computed and measured wall pressure distribution for the Bachalo and John- 
son transonic bump is given in Fig. 5. The incoming free stream is subsonic. A 
supersonic pocket is formed atop the bump as the flow accelerates through the 
bump. The combined effect of the shock-wave and the bump geometry results in 
flow separation near x/c = 0.7. KE2 modeP® shows the best overall agreement 
with the measurement than the three other models tested. KEl modeP^’^®, a low 
Reynolds number model, also gives better results than the Chien’s model. 

In Fig. 6, the temperature profiles at four measurement stations in the axial 
direction for the UC Irvine axisymmetric can combustor are shown. Compared with 
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previous published results^'^ obtained by using the KIVA code, the pdf chemistry 
model give a better agreement with measurement. The standard fc — e model was 
used in both calculations. 

2.2 Low-Reynolds Number Multiple-Sc€ile Model 

Damping functions were introduced into the high-Reynolds number multiple- 
scale model developed earlier^ so that the model transport equations can be inte- 
grated down to the wall. The Kolmogorov behavior of near-wall turbulence pro- 
posed by Shih and Lumley^^ were applied to determine the boundary conditions 
for the turbulent quantities at the wall. The details of the model can be found in 
a NASA TM in preparation. Fig. 7 shows the skin friction distribution for a flat 
plate turbulent boundary layer. The free stream Mach number is 2.87. The current 
model agrees better with the Van Driest II formula than the Launder and Sharma 
model. 

3. Future Plans 

(1) Examine the eifects of turbulence models on chemically reacting flow calcu- 
lations. Select a combustor of a different set-up, if necessary. 

(2) Broaden the validation cases to include realistic geometries found in gas tur- 
bine engine. This will provide means of realistically projecting model performance 
in engine component flows. 

(3) Continue the development of compressible turbulence models. 
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Figure 4. Compaxison of wall pressure for ONERA BUMP A. 



Figure 5. Comparison of wall pressure for the transonic shock/turbulent boundary-layer of 
Bachalo and Johnson. 
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Figure 6. Comparison of temperature (°K ) profiles for the UC Irvine axisymmetric can 
combustor. 
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Figure 7. Comparison of skin friction coefficients for a supersonic flat plate turbulent boundary- 
layer of Moo=2.87. 
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Modeling of Turbulent, Reacting 
Flows by PDF Methods 

A. T. Norris 

1. Motivation and Objective 

The objective of this work is the development, implementation and validation 
of Probability Density Function (PDF) models for turbulent, reacting flows. In 
addition, version 1.1 of the LPDF2D code was made available for release, containing 
many new features and able to be run in a UNIX environment. 


2. Work Performed 

In this section, the different areas of PDF methods that were worked on this year 
are described. Some vahdation tests axe also described. 

2.1 Introduction. 

The basis of Probability Density Function (PDF) methods hes in solving the 
transport equation for the joint PDF of quantities of interest, such as velocity, 
dissipation, enthalpy and/or composition. The reason why this approach to mod- 
ehng turbulent reacting flows is attractive is that the chemical reaction is treated 
exactly^ . 

Consider the transport equation for species (pa'- 


d(pa jj 

dt "dxi 


dJt 

dxi 


+ Sa\ 


( 1 ) 


where Ui is the velocity, Jf is the scalar flux and Sa is the chemical source term. 
Now the chemical source term is a unique function of the set of cr species and two 
state variables; 

Sa= Sa{(pl,-,4>a,P,T), ( 2 ) 

where p and T are the state variables pressure and temperature respectively. For 
all but the simplest of turbulent reactive flows, a solution of Eq.(l) is impossible. 
However the problem becomes tractable if the solution of the mean scalar field is 
attempted: 

+ + = + (3) 

dt dxi oxi oxi 

where (Q) denotes the mean of Q, and {Ui^a) is the velocity-scalar covariance, an 
unknown quantity that needs to be modeled. The other term that needs to be 
modeled is the mean reaction source term, (5 q). The reason that this needs to be 
modeled is because: 


Sai<Pu <Po,P,T), ^ Sa{{(pl), ( 0 .), {P), {T))). 


PRECEiUfU 



(4) 
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It has been shown^’^ that the two sides in this inequality can differ in both sign 
and by orders of magnitude. In addition, the modeling of the reaction source term 
would require different formulations for different fuels and mixtures. This would be 
clearly impractical. 

It is because of the difficulties in modeUng the reaction source term that PDF 
methods hold such attraction for the solving of turbulent reactive flows. 

2.2 Hybrid PDF Method. 

The PDF code under development is the LPDF2D solver developed by Hsu et 
al^. (1993). This code consists of a particle-based Monte Carlo solver for the 
solution of the transport equations for the PDF of the species and enthalpy, coupled 
with a finite-volume flow solver for the velocity, turbulence and pressure fields. 
Because this method is a combination of PDF and moment closure methods, it 
is referred to as a Hybrid PDF scheme. The solver can be used to solve steady, 
2D turbulent compressible reacting flows and has been extended to 3D. Work this 
year has focused on the refinement of the numerical algorithms, implementation of 
several different schemes to calculate reaction rates, porting the code to a UNIX 
workstation environment, parallel implementation of the code and some validation 
work. 

2.3 Numerical Refinement 

One of the drawbacks in implementing a PDF method is the extra computer 
memory required by the Monte Carlo solver. To help reduce the amount of memory 
required, important numerical developments were performed on the PDF solver: a 
new time- averaging was developed and implemented, and a convection algorithm 
that is accurate even with only a few particles was incorporated in the code. 

2.3.1 Time Averaging Scheme. 

One of the requirements of the Monte Carlo PDF solver is that it return smooth, 
or slowly varying values of the mean scalar to the finite-volume part of the code, 
otherwise the flow solver may not converge or could blow-up. The simple way to 
achieve this is to increase the number of samples being averaged over. However this 
will also increase the computer memory requirements. 

To achieve a mean scalar value with a small fluctuation, but with only a few 
particles, the technique of time averaging is used. The requirements of a good 
time averaging system are that only one set of data need to be stored, that the 
influence of samples from an earlier time can be m i n i m ized and that computation 
of the average is quick and simple. To fulfill these requirements, the time weighted 
average scheme was developed. In this method, the weighted time-average of some 
mean quantity {4>) at the nth time step is given by: 

(^)n ~ ~ (^) 

where the over-bar indicates a time averaged quantity and Wn is a weighting func- 
tion. 


Wn = Cn{Wn-l + 1 ). 


( 6 ) 
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where Cn is a function of n with a value between 0 and 1. For the case of c„ = 0, 
this scheme reduces to the case of no time averaging. With = 1, the average is 
formed over all time steps, each equally weighted. 

By trial and error it was found that a value of Cn = 0.9 provides a good compro- 
mise between quick response time and smoothing for the initial stages of a calcula- 
tion. When the solution of both the scalar and flow flelds has stopped exhibiting 
transient effects, an increase to c„ = 0.999 reduces the scalar fluctuation further, 
allowing the solution to converge. 

It should be noted here, that the concept of convergence in PDF methods is not 
the same as that in pure moment closure schemes. Due to the statistical error 
always present, solutions can only converge down to the level fluctuation present in 
the mean estimates, rather than to machine ciccuracy. 

2.3.2 Convection Scheme. 

Convection in the Monte Carlo PDF code consists of the movement of particles 
from one cell node to another due to the action of the mean velocity and turbulent 
diffusion. Due to the finite number of particles at each node and the need to move 
whole numbers of particles, errors occur in this process. For example, if 5.3 particles 
are to be convected, the existing code would convect 5 particles. By increasing the 
number of particles at a node, more will be convected and the errors would be 
decreased, however more particles means more memory required. 

To overcome this problem, we make use of the time-averaging used in the code. 
At each step, the number of particles to be convected is calculated. This number 
is then rounded to an integer value by a random process, with the probability of 
rounding up or down dependent on the value of the fractional remainder. ^From the 
previous example, 5.3 has a fractional remainder of 0.3. Thus 5.3 would be rounded 
up to 6 with a 30% probability, and rounded down to 5 with a 70% probability. 
Because of the time averaging used in the scheme, the mean convection will be 5.3 
particles. 

Another method would be to save the fractional particles and carry them over to 
the next time step. This however has two disadvantages. First the fractional value 
needs to be saved, thus using memory. Secondly, this introduces an occilation into 
the solution. For example, if 5.1 particles were to be convected, then 5 particles 
would be convected except for every 10th step when 6 particles would move. This 
period would be damped out by the time averaging to some extent, but not totally. 
This small periodic fluctuation could possibly excite some non-physical fluctuations 
in the solution, such as vortex shedding. For these reasons, this method was not 
implemented. 

2.4 Operating Environments. 

In order to make the LPDF2D code accessible to a wider range of users, several 
changes have been made to the code. The most common computing environment 
available these days is the UNIX workstation. Because of this, the LPDF2D code 
has been modified to run on workstations. This has involved the removal of machine 
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dependent commands, the introduction of a portable random number generator and 
the ability to perform calculations in a parallel cluster environment. 

2.4 -1 Random Number Generator. 

To aid the portability of the code, a machine independent random number gener- 
ator was incorporated into the code. This routine provides a "more random” set of 
numbers than most built-in generators, while making use of the 32 bit wrap-around 
feature of UNIX machines to speed up the program. 

2.5 Chemical Kinetics. 

Despite the ability of the PDF method to treat chemical reaction exactly, the 
implementation of the numerical chemical kinetics in a Monte Carlo scheme is a 
non-trivial matter. First, the solution of the full system of rate equations for the 
thermochemistry in turbulent reacting flows requires obtaining solutions for of or- 
der 50 chemical species, governed by of order 200 stiff, non-linear rate equations. 
At present this task is computationally infeasible, and so reduced mechanisms are 
employed. These reduced mechanisms represent the full composition by a few rep- 
resentative species, typically two to five, and the reaction rates of these species are 
quasi-global equations derived from the full mechanism. However even the solution 
of a reduced mechanism is a computationally expensive task. In a typical PDF 
calculation®, there are on order 100,000 particles and 1,000 time steps, resulting in 
the code performing 10^ solutions of the reduced mechanism. 

A computationally eflicient way of reducing this task to manageable proportions 
is to use a look-up table. In this table, the reduced mechanism is integrated for 
discrete time and composition increments, and stored in a table. Thus integration 
is replaced by interpolation, at a considerable saving in computer time. By the use 
of adaptive tabulation techniques'^ the size of the tables can also be minimized. 

However, because of the many different combinations of fuel and oxidizer, as well 
as the differing degrees of complexity of mechanisms, several different options for 
chemical reaction have been provided apart from look-up tables. These include 
equilibrium chemistry (where the reaction is assumed to proceed infinitely fast), 1 
step global reactions and the option to integrate the full mechanism. 

In these options, the CHEMKIN® package has been used to provide the species 
data for the reactions, reducing the amount of input required by the user. 

2.6 Parallel Computations. 

If the option to use the full mechanism is used, the biggest use of time in the 
PDF code is the calculation of reaction rates. Because of this, the LPDF2D code 
has been altered to run on a work-station cluster, with PVM'^ message passing. At 
present this involves only the calculation of reaction rates in parallel, the rest of the 
process being performed in serial. Because of the low volume of message passing, 
and the statistical nature of the PDF code, this has proved to be a very effective 
method of achieving substantial speedup of the code. 

4.0 Validation. 

In order to validate the PDF code, and test its performance against that of a 
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conventional moment closure code with laminar chemistry, a test case was chosen 
and then modeled by both methods. The flow chosen was the CO/H 2 /N 2 -air flame 
of Masri et al®. This flow was chosen as it has high levels of turbulence in the 
reactive region of the flame, thus the turbulent /chemical interactions would have a 
significant effect on the flow. 

The conclusion of the study was that the PDF method did provide a superior 
performance to the moment closure method in predicting turbulent reacting flows. 
In Figs. 1 and 2, the temperature contour plots for the two numerical methods are 
shown. It can be seen that the PDF method predicts peak temperatures of about 
1400K in a narrow band, in good agreement with the experimental data. However 
the moment closure results show peak temperatures of 2000K lying in a broad band, 
contrary to the experimental results. 

Details of this comparison are described by Norris and Hsu^, and were presented 
at the 30th Joint Propulsion Conference in Indianapohs. 
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Figure I: ^ntour plots of mean flame temperature produced by the PDF method. Note that radial 
coordmate is stretched by factor of five. Contour levels are in degrees Kelvin. 
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Application of PDF to Compressible 
Turbulent Reactive Flows 

Andrew T. Hsu 


1. Motivations and Objectives 

In an earlier paper[l], we have introduced a probability density function (pdf) 
approach for high speed reactive flows. The extension of this pdf approach to three- 
dimensional applications is reported here. The major obstacle to the application of 
a Monte Carlo pdf solver to complex three-dimensional reactive flows is the limita- 
tion imposed by the need of a large number of notational particles in a Monte Carlo 
simtilation to provide a converged, relatively smooth solution free of large statistical 
errors, which in general requires enormous amount of both computer memory eind 
cpu time. To overcome this obstacle, a novel averaging scheme is introduced. This 
new averaging scheme allows one to use as few as 50 sample particles per compu- 
tational cell and yet still provides Monte Carlo solutions that are smooth enough 
to be applicable in the finite-volume/ Monte- Carlo scheme. Further, the averaging 
scheme automatically eliminates the effect of initial conditions, mak ing it ideal for 
large-scale applications. 

2. Work Accomplished 

In turbulent reactive flow computations, the conventional moment closure models 
have difficulty in treating the nonlinear chemical reaction source terms; this problem 
is known as the chemical closure problem. One of the methods that csin overcome the 
chemical closure problem is the probability density function (pdf) method. Several 
key advances were made in the development of pdf methods during the past decade. 
Most of the earlier advances are primarily in the low-speed combustion area. 

In an earlier paper a probability density function turbulence model for compress- 
ible reacting flows has been proposed by the present authors [1]. The probability 
density function of the species mass fraction and enthalpy is obtained by solving a 
pdf evolution equation using a Monte-Carlo scheme. The pdf solution procedure is 
coupled with a compressible finite volume flow solver which provides the velocity 
and pressure fields. A modeled pdf equation for compressible flows, capable of cap- 
turing shock waves and suitable to the present coupling scheme, has been proposed 
and tested. Several 2D supersonic diffusion flames were studied and the results 
compared favorably with the available experimental data. 

In the present paper, the extension of the pdf method to three-dimensional su- 
personic combustion is reported. The main obstacle to the application of a Monte 
Carlo pdf solver to complex three-dimensional reactive flows is the limitation of 
computer memory and cpu time. To overcome this obstacle, a novel averaging 
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scheme is introduced. This new averaging scheme allows one with the use as few as 
50 sample paxticles per computational cell to obtain Monte Carlo solutions that are 
smooth enough to be supplied to the finite-volume solver as required during each 
of the iteration steps of the finite- volume/Monte-Carlo scheme. Furthermore, the 
averaging scheme automatically eliminates the effect of initial conditions, maldng 
it ideal for large-scale steady-state applications. 

A three-dimensional jet-in-crossfiow hydrogen combustion case has been success- 
fully computed using the present scheme. The results are compared with those from 
a finite volume reactive flow solver. 

2.1 The PDF Method 

The modeled pdf equation for the species mass fraction and specific enthalpy used 
in the present work can be written as follows: 

{pP),t + (/>< > P)j + (pt->jP),Y, = (DtP.j)a + M(P) - (SpP),». 

where the first two terms on the right hand side of the equation are the mod- 
eled terms for turbulent diffusion and molecular mixing; the last term is the term 
representing the compressibility effect. 

The turbulent diffusion term is simulated using a simple gradient diffusion model 
by 0’Brien[2] and Pope [3]: 


-<u)\Yi,h>P = DtP^j. 


The use of this gradient diffusion model is consistent with the use of a A: — e model 
in the flow solver. 

The molecular diffusion term is modeled using the continuous mixing model by 
Hsu and Chen [4]. This model is an extension of Curl’s model. In order to achieve 
continuous mixing, we a.ssume that all the particles within a cell participate in 
mixing. The extent of the mixing is controlled at the individual particle level. That 
is to say, the N particles within a given cell are randomly grouped into N /2 pairs; 
the properties of all the particles change according to 


+ et) = AYm{i) + {1- A)Yn(t) 

Ym{t + St) = AYr^it) + (1 - A)y^(t) 

The extent of mixing is controlled at the individual particle level through the pa- 
rameter A, which is defined as 

, St 
A = C'^- 

T 

where C = 2.0. At the hmit St — > 0, the above equations becomes 

T 


dt 
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The above equation states that the change of Yn due to mixing is proportional to 
the difference between Ym and Yn, and inversely proportional to the turbulence time 
scale, T. 

The compressibility effect is modeled by the following expression: 

Sp =< p>,t + < Ui X p >,i +0.8/9 < k >< Ui >,,• —a2pPrMt + azpeMf , 

which caji be regarded as the convection velocity of a sample particle in the h- 
direction in the space spanned by h and Ti’s. Details of this model are given in Ref. 
[ 1 ] 

2.2 Solution Procedure 

A fractional- step Monte-Carlo method developed by Pope [3] is used in solving 
the pdf evolution equation. When coupling a pdf solver with RPLUS, a finite 
volume flow solver, the species transport equations in the RPLUS code are no 
longer needed. The information we need, at each marching time step, from the flow 
solver (RPLUS) includes the mean velocity, pressure, density, and a turbulence time 
scale or quantities such as the turbulent kinetic energy and dissipation rate. The 
species transport and chemical reactions are simulated by solving the pdf evolution 
equation. At every time step, the temperature calculated from the pdf solution is 
fed back to the mean flow solver for the computation of heat transfer and pressure. 
The Monte Carlo and finite-volume solvers are run in parallel, and information 
exchange occurs at every time step until a converged solution is obtained. 


2.3 Averaging Scheme 

In a Monte Carlo computation, the statistical error is estimated to be of the 
order of I/VN, where N is the number of sample particles per computational cell. 
If this error is large, it could lead to an erroneous finite-volume solution when 
information from the pdf solver is transferred to the flow solver. Because of this 
slow convergence, to obtain a smooth solution from the Monte Carlo solver, one 
needs thousands or even tens of thousands of sample particles per cell. In our 
previous paper [1] we proposed a combined ensemble-time averaging scheme. That 
scheme reduced the required particle per cell, but it requires some more memory 
for the time averaging process. For any realistic 3D combustion problems, the 
memory requirement of that scheme exceeds the CRAY-YMP capacity. In what 
follows we introduce a novel averaging scheme that eliminates the effects of the 
initial conditions automatically while requiring no additional memory. 

Let a„ be the combined time-ensemble average at the time step, and let bn 
be the ensemble average at the time step. A general weighted time average can 
be written as ^ 

a^n — "7 

Jn t J- 

In terms of ensemble average, the above equation can be written as 


an = 


fn 


■hn-l 


f-n.f-n.-Z 


/. + 1 (/„ + + 1) (/„ + l)(/..-l+l)(/„-2+l) 


^n-2 + ‘ 
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Where/, rs e discrete function that can be chosen appropriately to give the desired 
weighting for so utions from various time levels. For instance, /, = n would give a 
equaUy weighted miming average, and the above equation becomes 


= 


n + 1 


(&n + hn-i + f- + 6o) 


Fig 1 shows the weighting distribution as a function of the time step when using 
tne following function ® 

/n = nc 

jrith c = 0 6, 0.8, 10. This weighting procedure gives a fairly smooth Monte 

^^osing^"" ^ ehminates the effect of only the first few time steps. By 

/n ~ ^{fn—1 + 1) 

In which case Eq. (9) is equivalent to 


ar, = 


+ + • • • + c^-% + c"6o) 


The weighting distribution, as shown in Fig. 2, is shifted drastically towards the 

latest solutions The weakness of this particular formulation is unevenness of the 
weight distribution. 

IdeaUy, one wo^d like to eliminate the effect of the initial conditions, yet have 
a fairly even weight distribution for the later solutions. Such a weighting function 
can be achieved by replacing the constant c in the equation for /„ by a variable c„. 


fn — ^n(/n-l + 1) 


Th«e are any number of possibilities in choosing c„. For example, c,, = 1 - 1/n 

r”. , ~ ^ ~ formulation we have found to be 

leurly successful is 

Ti!*' j ^ V distribution for x = 1.014, 1.0088, and 1 007 

roud^’p of this waghted averaging procedure are that the effect of the initial 

condition can be effectively ehminated, and that after a certain number of iterations 
the solutions are equaUy weighted. For example, in the case of x = 1.014, the effect 

lout Sott 'r t and after 

about 500 time steps, the solutions are equally weighted in the averaging procedures. 

The pdf results reported in what foUows are obtained using this weighting scheme. 
2.4 3D Application 

In the present work, one of our major objectives is to demonstrate the feasibility 
of applying the pdf method to realistic 3D supersonic combustion problems. Unfor^ 
tunately, there is a lack of experimental data in this area. The case we have chosen 
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to demonstrate the application of the pdf method is a supersonic jet in crossflow. 
A 300°iiir, Mach 1.5 hydrogen jet is vertically injected into a Mach 2 crossflow with 
1000°ir vitiated air. The jet diameter is 6 mm. Figs. 4-8 are the temperature 
and species mass fraction contours from the pdf solution as comp£Lred to the solu- 
tion from a finite volume code (RPLUS) using laminar chemistry. The grids used 
in both calculations are 40 x 15 x 25. The computational domain reaches about 
30 jet diameters downstream of the jet. The Monte Carlo simulation of the pdf 
equation used 50 sample particles per computational ceU. The coupling between 
the pdf solver and the finite volume flow solver is done through feeding density, 
velocity, and turbulent time scale from the finite volume solution to the pdf solver, 
and feeding temperature distribution from the pdf solver to the finite-volume flow 
solver. 

Fig. 4 shows the central-plane temperature distribution. The results show that 
the pdf solution gives a much shorter cold jet core than that from the finite volume 
scheme, which is consistent with our previous 2D results. The majdmum tempera- 
ture from both solutions are the same. Fig. 5 shows the temperature contours at 
a cross section at about 12 jet diameters downstream of the injection port. The 
basic structure of the two solutions axe similar, with some differences in details: 
The finite volume solution gives a very high temperature region at the bottom wall 
away from the centerline while the pdf solution does not show this. This difference 
could be caused by the different boundary conditions used at the wall. 

Figs. 6, 7, and 8 show the mass-fraction contours of hydrogen, oxygen, axid water 
vapor, respectively. The results for the mass fractions are consistent with what 
we have observed from the temperature distribution. For instance, the hydrogen 
distribution from the finite volume solution shows a maximum of about 0.8 near the 
center while the pdf solution there has a lower value of only 0.35 indicating some of 
the differences observed between the solutions. In Fig. 8, where it corresponds with 
the high temperature at the wall, there is a higher concentration of water vapor 
(about 0.35) in the finite- volume solution while the pdf solution shows a maximum 
of only 0.12. 

Because of the lack of experimental data, it is difficult to judge the validity of 
either solution. However, the contours from the Monte Carlo solution is as smooth 
as that from the finite volume solution, showing that it is fecisible to malce large 
scale computations on practical problems using the Monte Carlo pdf solver. 

2.5 Concluding Remarks 

A pdf method for high speed turbtilent combustion hats been successfully extended 
to 3D applications. A novel averaging scheme has been introduced which made the 
applications of the present pdf method to realistic large scale computations feasible. 

3. Future Work 

3.1 PDF 

Parallel computing and NOx prediction will be the two major areas for future 
PDF applications. The goal is to make the PDF method an industrial design tool. 
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Fig. 2. Weight distribution for the first 1000 time 
steps using Eq. (12). 



Fig. 1. Weight distribution for the first 1000 time 
steps using Eq. (11). 



Fig. 3. Weight distribution for the first 1000 time 
steps using Eqs. (14) and (15). 
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Fig. 4. Center plane temperature distribution for the 
3D supersonic jet combustion problem. 
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Fig. 5._ Temperature distribution at x=12d for the 
3D supersonic jet combustion problem. 
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Turbulence Model Development for NPARC 

J. Zhu, Z. Yang, and T. H. Shih 

1. Motivation and Objective 

The purpose of this research is to enhance the turbulence modehng capability of 
NPARC. The outcome of this research is a turbulence subprogram which is indepen- 
dent of NPARC but linked to NPARC in code execution. In the solution process, 
the turbulence subprogram interacts with the NPARC code, the mean field solver, 
to produce the solution. This turbulence subprogram contains a number of turbu- 
lence models suitable for aerospace and aero-propulsion system appUcations. These 
models axe chosen from the state of the art of turbulence models developed by both 
CMOTT researchers and researchers from other groups of the turbulence modeling 
community. The turbulence subprogram with advanced turbulence models serves 
the following purposes; first, it will enhance the turbulence modeling capability of 
NPARC; second it will provide a vehicle for incorporating future turbulence mod- 
els, and third it can be used to validate turbulence models against complex flows of 
industry interest. 

NPARC is a general purpose computational fluid dynamics (CFD) code. It solves 
Navier-Stokes equations in a conservation form in a general coordinate system. It 
can handle complex geometries and different types of boundary conditions. The 
NPARC code is widely used in the aerospace and aero-propulsion community. Up 
to now, many flows of practical interest have been calculated using the NPARC 
code, ranging from an inlet flow in an air-breathing engine to the flows around a 
maneuvering airplane. A description of the code and its numerical scheme is given 
by Cooper and SirbaughU Updates of code capabihties and a Usting of bibhogra- 
phies of using the NPARC code to calculate flows of practical interests can be found 
in the current version of the NPARC User Guide^. The NPARC code is actively 
supported by NASA Lewis Research Center via the NPARC Alliance, a partnership 
formed between NASA Lewis Research Center and Air Force Arnold Engineering 
Development Center. The NPARC users are represented by the NPARC Association 
currently headed by Boeing Company and McDonnell Douglas Company. 

2. Work Accomplished 

To enhance the turbulence modehng capabihties of NPARC, a turbulence module 
approach is adopted in the present project. In the foUowing, the turbulence module 
for NPARC is briefly described. We wiU then hst the turbulence models that are 
available in the current version of the turbulence module. NPARC calculations 
using the turbulence module for some propulsion flows will also be presented. 

The work accomphshed is done by J. Zhu, Z. Yang, T. H. Shih at Center for 
Modehng of Turbulence and Transition (CMOTT), and by N. Georgiadis at NASA 
Lewis Research Center. This project is coordinated by D. R. Reddy of NASA Lewis 
Research Center. 
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2.1 Introduction to Turbulence Models in NPARC 

Most of the flows in aerospace and aero- propulsion systems are turbulent because 
of the high Reynolds numbers of the flows. Thus, for a code with a well developed 
numerical scheme such as the NPARC code, its turbulence modeling capability is 
the pacing item for computational accuracy. Currently, NPARC contains a num- 
ber of mixing length eddy viscosity models^’^ and the Chien k — e eddy viscosity 
modeP. A mixing length eddy viscosity model has the advantage of being simple in 
implementation and fast in execution; however it is a ‘one physics/one flow’ model, 
i.e., the mode is tailered to a particular physics/flow, and its performance deterio- 
rates drastically when used to calculate other flows. In comparison, & k — e eddy 
viscosity model is much more general than mixing length eddy viscosity models 
and could be used for a much wider range of flows, for example, both free shear 
flows and wall-bounded flows. Among the existing low-Reynolds number k-e eddy 
viscosity models, the Chien k — e eddy viscosity model is perhaps the simplest and 
gives a reasonable performance for attached boundary layer flows. However, it has 
some well known deficiencies. For example, its performance is rather poor for flows 
with transition, separation and reattachment, because is used in modeling the 
near- wall effect. These deficiencies caji be removed and the applicability of A: — e 
model can be substantially expanded. Recently, a CMOTT k-e eddy viscosity 
model has been implemented into the NPARC code. This model does not contain 
the friction velocity, unlike the Chien k-e model, and can thus be used for flows 
with separation and reattachment. The detail of the model form and the results of 
model performance will be presented in later subsections. 

2.2 CMOTT Approach - Turbulence Subprogram for NPARC 

In the effort of enhancing the turbulence modeling capabihties of NPARC, we 
implement different turbulence models in a stand-alone turbulence subprogram (also 
referred to as turbulence module) instead of implementing them in the NPARC code 
directly. We feel that this is a preferred approach in a situation that a gap exists 
between turbulence model developers and CFD users. The former mainly use simple 
flows to verify new modeling concepts and evaluate the resulting models, while the 
latter axe usually reluctant to implement new turbulence models unless they see 
that such models have shown a good performance for a wide range of complex flows 
relevant to industry situations. In order to bridge this gap, we develop turbulence 
modules for industry’s CFD codes. Since the NPARC code is a major industry code 
with a well established numerical scheme and a large user group in the aerospace 
and aero-propulsion community, we use the NPARC code as our reference CFD 
code to write the turbulence modules. With minor modifications, the turbulence 
module written for the NPARC code can be adopted for other industry CFD codes. 

The turbulence module is written in a self-contained manner so that the user can 
use any turbulence model in the module without worrying about how it is imple- 
mented and solved. The inputs to the turbulence module are mean flow variables, 
boundary and geometric information which are provided by a mean flow solver. The 
outputs of the turbulence module are the Reynolds stresses and heat fluxes which 
are needed for the mean flow calculation. In the current formulation, the Reynolds 
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stresses and heat fluxes are expressed in terms of turbulent diffusivities and rel- 
evant turbulent source terms. The interaction between the mean flow solver and 
the turbulence module gives the final turbulent flow solution. With the aid of the 
turbulence module, we can take the advantage of the well-established and sophisti- 
cated CFD codes to test turbulence models for a variety of complex flows which are 
intractable with the simple research codes. The turbulence module, on the other 
hand, enhances the turbulence modeling capabilities of existing CFD codes and 
provides a vehicle to facilitate the technology transfer from the model development 
to model applications. 

In the turbulence subprogram for NPARC, the governing transport equations 
for turbulence quantities {k and e) are written in conservation form and a finite- 
volume procedure is used to solve these equations. Comparing with the NPARC 
code, the turbulence subprogram has the following features; (a) the cell-corner 
arrangement is used, i.e., the flow variables are stored at grid nodes rather than at 
the centers of control volumes; (b) contrary to the NPARC code and most other 
codes for compressible flows, the non-delta form of equations is used which leads 
to a simpler linearization and is more effective to ensure the positiveness of the 
turbulent kinetic energy and its dissipation rate (fc and e) than the delta form. 
To reduce numerical diffusion while maintaining necessary stability, the second- 
order accurate HLPA scheme® is used for the convective terms of the turbulent 
transport equations. This scheme is implemented such as to blend two schemes by 
a parameter A, with Umiting values A=0 for the (first-order) upwind and A=1 for the 
second-order and bounded solution scheme. The turbulence subprogram includes 
four most commonly used boundary conditions - the infiow condition, the outflow 
condition, the symmetry plane condition, and the condition on the solid waU. The 
turbulence subprogram also allows the user to introduce boundary conditions other 
than these four types. Care has been taken to make the turbulence module user- 
friendly. This is achieved by providing a few parameters that the user can define in 
model selection and problem specification. The resulting algebraic set of equations 
are solved iteratively via a under-relaxation procedure. The user can also control 
the relaxation parameters in the problem specific part of the subprogram. The 
turbulence module is written in a similar manner as the FAST2D code written by 
Zhu^. A detailed description of the turbulence module for the NPARC code is 
going to be presented in a paper by Zhu and Shih® in the upcoming 31st AIAA 
Joint Propulsion Conference. 

2.3 Turbulence Models in the Subprogram 

Currently, four turbulence models are available in the turbulence subprogram. 
They are the mixing length eddy viscosity model of Baldwin-Lomax®, the low 
Reynolds number k - e eddy viscosity model of Chien®, the low Reynolds num- 
ber k — e eddy viscosity model of Shih and Lumley^, and a new low Reynolds 
number reahzable k - e eddy viscosity model developed here at CMOTT. 

Because of the increasing computer power available and the increasing demand 
for more accurate flow prediction, it is recommended that a A; - e eddy viscosity 
model is used for calculating flows of engineering interests. The Baldwin-Lomax 
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mixing length eddy viscosity model is included in the turbulence subprogram for 
two reasons. First, it is chosen as a representative of the mixing length models. 
Second, it can be used to generate the initial field for k — e model calculations. 

The Baldwin-Lomax mixing-length model, the Chien k-e model, and the Shih- 
Lumley k-e model have been published in archival journals. Here, only the newly 
developed low Reynolds number realizable k-e eddy viscosity model is described. 
The model has the same model form as the Shih-Lumley low Reynolds number 
k-e eddy viscosity model, i.e, it has the same form of transport equations for 
k and e, it uses the same near-wall boundary conditions for k and e as derived 
from the Kolmogorov behavior of the turbulence field near the wall, and it uses 
the same parameter Ry = yk^/^ /v to construct the damping function /^. The 
boundary conditions based on Kolmogorov behaviors of near wall turbulence makes 
the computing robust. The parameter Ry is a better choice for modeling of near 
wall turbulence than ?/+ since it does not involve the skin friction Ut- The major 
improvement in the present realizable k-e eddy viscosity model compared with the 
Shih-Lumley model is that a realizable formulation for is used. In an isotropic 
eddy viscosity model, the Reynolds stresses are related to the mean velocity field 
by 

2 2 

UiUj -)- Uj^i — —Uk,k^ij') ( 1 ) 

In the framework of fc — e two equation eddy viscosity model, the eddy viscosity i/y 
is further expressed as 

I'T = C'm/m Y (2) 


where the damping function is introduced to account for the effects of wall. 

Traditionally, a constant value of = 0.09 is used in equation (2). However, 
any two-equation eddy viscosity model with a constant is unrealizable in the 
sense that the turbulent component energy given by the model can become nega- 
tive in certain flow situations. In addition, data from direct numerical simulations 
and from experiments suggest that vary with fiow situations. For example, to 
match the data, different values should be used for homogeneous shear fiows 
and for turbulent flows in a two dimensional channel. In this study, a variable 
formulation is used, in which depends on the solutions of the mean flow field 
and the turbulence field. The formulation was derived based on the realizability 
constrains, rapid distortion theory (RDT), and the invariants theory and takes the 
following form 




1 

Ao + A,U*!i' 


( 3 ) 


In equation (3), k/e represents the characteristic time scale associated with the 
turbulent motion. U* is related to the characteristic time scale of the mean field 
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and it is given by 

V = 

Stj = Sii - 

fijj- = D,^j — 2eijkU>k 

where Sij is the mean strain tensor, Q,*j is the mean rotation tensor in the inertia 
frame, and Uk is the angular velocity of the rotating reference frame. The mean 
rotation tensors in the inertia frame and in the rotating frame are related by 


fljj — fijji ^ijk^k- (^) 

In a rotating frame, the mean strain tensor Sij and the mean rotation tensor Clij 
are related to the mean velocity field by 

Sij = l{Ui,j + Uj,i), Qij = ^{Ui,j - Uj,i). ( 6 ) 

It is worth noting that the effect of frame rotation is included in the present model 
directly via its effect on U*, and consequently on and 

The parameter ^4^ in equation (3) is a scalar function formed from the invariants 
of the mean velocity field. It has the following expression 


As = VQcos(f>, (p = -arccos{V6W) 

O 


( 7 ) 


As shown in Reynolds^°, the choice of A^ guarantees the correct hmiting behavior 
of the model in the rapid distortion limit. 

The parameter Aq in equation (3) is a model constant. In the present study, 
Aq = 4.0 is chosen so that the present formulation gives a value of 0.09 
for the equihbrium log layer regions of two-dimensional turbulent boundary layer 
flows. It was based on the experimental data for these flows that a constant value 
of = 0.09 was first suggested. 

Equations (3) — (7) give the formulation that we use in the present study. 
For such a formulation, the resulting Reynolds stresses are always realizable in 
that no physically un-admissible Reynolds stresses will be produced. In addition, 
the variations of with flow situations are in the right direction, i.e., the value of 
Cfj. is smaller than 0.09 in regions where it should be. Thus, it can be expected that 
the current formulation of will improve the performance of k - e eddy viscosity 
model considerably, particularly in regions where strain rate is large. 

In a paper by Shih et al.^\ the above mentioned formulation was used in 
conjunction with a new (high- Reynolds number) k — e eddy viscosity model, in 
which the dissipation rate equation was derived based on the asymptotic behav- 
ior of the fluctuating vorticity equation in the hmit of large Reynolds numbers. 
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The paper shows that the Cfj, formulation helps the model performances for rotat- 
ing homogeneous shear flows, free shear flows, and wall boundary flows. For wall 
bounded flows, the near-wall boundary conditions were provided by the standard 
wall function formulation presented in Launder and Spalding^ 

2.4 Flows Calculations using NPARC and Turbulence Module 

A series of flow calculations have been made using the NPARC code to solve the 
mean flow field and the turbulence module to solve the turbulence model equations. 
The turbulence module was first validated against the turbulent boundary layer 
flow over a flat plate with zero pressure gradient. Complex flows of aero- propulsion 
interest were then calculated. These include an ejector nozzle flow which provides a 
good test for flow mixing; a transonic diffuser flow which tests shock wave/turbulent 
boundary layer interactions; and a boat-tail nozzle flow which tests flow separation 
and jet-boundary layer interactions. In each case, the numerical convergence and 
the grid sensitivity are checked to make sure that the numerical solutions obtained 
are creditable. In the following, each flow case is briefly reported. These calculations 
are going to be presented in a paper at the upcoming 31st AIAA Joint Propulsion 
Conference by Yang et al.^^. 

2. 4-1 Flow over a flat plate 

The purpose of calculating this flow is for code validation. The computational 
domain is shown in figure 1. The total grid points are 111 by 81 with 111 grid 
points in the x (streamwise) direction and 81 grid points in the y direction. In the 
X direction, 14 grid points are located before the leading edge of the flat plate. In 
the y direction, grid points are stretched so that the first off wall points have an 
average value of about 1. Since the NPARC code is for compressible flows, and 
computational results are going to be compared with an incompressible turbulent 
boundary layer, a freestream Mach number of 0.2 is used in the present calculation. 

Although straightforward on the surface, the calculation of this flow is quite 
demanding in that the flow starts from a laminar-like flow, then undergoes a tran- 
sition process, and then finally settles down to flow of the turbulent boundary layer 
type. It is very difficult to capture all three flow regions accurately. In the present 
computation, we are mainly interested in the solutions in the turbulent boundary 
layer region. Figure 2 shows the skin friction coefficient distribution as a function 
of Ree, the Reynolds number based on the momentum thickness of the boundary 
layer. Both the model predictions and the experimental results of Wieghardt and 
Willmann^^ are shown for comparison. It is seen that all models give similar results 
and all model predictions are close to the experimental data. Comparisons were 
also made between solutions of the NPARC code with built-in turbulence models 
and solutions of the NPARC code with the turbulence module. These solutions are 
found to be the same, thus the code with the turbulence module is validated. 

2-4 -2 Flow in an ejector nozzle 

Ejector nozzle flows have important applications in propulsion systems because 
ejector nozzle is a candidate device to reduce the jet noise from the engine. The 
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noise reduction is particularly relevant to the High Speed Research (HSR) Program 
and the High Speed Civil Transport (HSCT) Program in NASA. A decade or two 
ago, ejector nozzle was also studied for its applications in the vertical land/taking 
off aircraft. 

The geometry for an ejector nozzle is shown in figure 3. The performance of 
the ejector nozzle is measured by its mixing properties. To see the mixing effect, 
computed velocity and temperature profiles at a few stations along the streamwise 
direction are shown in figure 4 and figure 5, respectively, along with the experimental 
data of Gilbert and HilP^ It is seen that different k-e eddy viscosity models give 
similar results and that results from k-e eddy viscosity models are much better 
than those from the mixing length eddy viscosity models. The same conclusion was 
also reached by Georgiadis and Yoder^® who calculated this flow using the NPARC 
code with its built-in turbulence models. 

2 . 4.3 Flow in transonic diffuser 

Transonic diffuser flow is another flow of interest in propulsion systems. In the 
present calculation, the experiment by Salmon, Bogar, and Sajben^^ is chosen as 
the test case. Depending on the pressure at the nozzle exit, there can be three 
different kinds of flow patterns. They are the no-shock case, the weak-shock case, 
and the strong-shock case. All these three cases were calculated in the present 
study, however, only the results for the strong-shock case, which is also the most 
difficult case, are shown here. 

Figure 6 shows the diffuser geometry, together with grid information. The pres- 
sure distribution along the upper wall of the diffuser is shown in figure 7. Both the 
experimental results and the results from different model predictions are shown. 
For this flow, the critical parameters that are of interest to the design engineers 
are the location of the shock wave and the pressure recovery after the shock wave. 
Among the model predictions, the new realizable k-e eddy viscosity model gives 
the best performance in these critical parameters. Figure 8 shows the velocity pro- 
file at A = 2.88, a station located after the shock wave where the experimental 
data are available. The realizable k-e eddy viscosity model gives a much better 
prediction for both the overall feature of the field and the size of the separation. 
Since the controlhng mechanism for this flow is the shock wave/turbulent bound- 
ary layer interaction, we can conclude that the new reahzable k-e eddy viscosity 
model captures this mechanism much better than the k-e eddy viscosity model 
with constant C^, such as the Chien k-e eddy viscosity model. The Chien k-e 
model is itself much better than the Baldwin-Lomax mixing length eddy viscosity 
model for this flow. 

2.4-4 Flow around a boat-tail nozzle 

Boat-tail nozzle is a device commonly found in propulsion systems. A particular 
interest of the present study is on the effect of the jet plume on the turbulent 
boundary layer formed on the nozzle afterbody. Depending on the flow situations, 
the boundary layer can be either separated or attached. Boat-tail nozzle flow also 
provides a good test case for the two dimensional version of the NPARC code in 
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that it tests the code’s capability for axisymmetric flows. 

The experimental studies of Mason and Putnam^*, Abeyounis and Putnam^® are 
chosen to validate the model calculation. In these studies, the freestream Mach 
number is 0.8. The stagnation pressure in the nozzle is such that the jet exits the 
nozzle at the local sonic speed, i.e. at Mach number of one. The nozzle geometry 
and the computational domain is shown in figure 9. Figure 10 shows the pressure 
distribution on the nozzle afterbody from both the computation and the experiment. 
It is seen that both the Chien k — e model and the new realizable eddy viscosity 
model give similar results, and they are all fairly close to the experimental data. 
The k—e eddy viscosity models give improved predictions compared with the mixing 
length eddy viscosity model. 

3. Future Plans 

Turbulence model development for NPARC is a long term activity. What was 
reported above represents the first part of this project. Many other activities need 
to be carried out in this project. For example, the numerical capability of the 
turbulence module needs to be enhanced; more axivanced turbulence models need 
to be incorporated into the turbulence subprogram; and NPARC calculations with 
the turbulence modules need to be carried out for a large range of flows relevant 
to propulsion systems. To this end, a research team has been formed. In addition 
to the CMOTT members for this project (J. Zhu, Z. Yang, and T.H. Shih), this 
research team also includes N. Georgiadis and D.R. Reddy of NASA Lewis Research 
Center, and J.R. Sirbaugh at NYMA Inc. The activities of this research team has 
been incorporated into the NPARC Development Plan of the NPARC Alliance. The 
following are some of the planned highlights of these activities. 

3.1 Turbulence Module for 3D Flows 

Currently, the turbulence module is written for 2D flows (both planar flows and 
axisymmetric flows.) Two dimensional flows provide a basis for code and model 
validation. However, a code for 2D flows is only the first step in the coding process 
since almost all practical flows in aerospace and aero-propulsion systems are three 
dimensional. Thus, the next step is to extend the turbulence module to three 
dimensional flows and to test the turbulence module and models accordingly. It is 
expected that a 3D version of the turbulence module will be available by Summer 
1995. 

3.2 Tensorial and Galilean Invariant Eddy Viscosity Model 

So far, all the k — e eddy viscosity models in the turbulence module depend 
on geometry information in their near wall treatment, e.g. the distance from wall 
appears in the damping function. Models with geometry information are undesirable 
for flows with very complex geometries because in those situations, the definition of 
wall distance can be ambiguous. In addition, models without geometry information 
can be easily incorporated into CFD codes using unstructured grids while models 
with geometry information can not. To propose such a tensorial invariant model, it 
is necessary to use parameters other than Ry or y+ to calibrate the near wall effect. 
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In the paper by Jones and Launder^*^, in which the first low-Reynolds number k-e 
eddy viscosity model was proposed, Rt, the turbulence Reynolds number, was used. 
In a more recent paper by Yang and Shih^^, a new parameter, R = k/{Su), was 
introduced to model the near wall effects. Both Rt and R are field variables and are 
candidates to construct tensorial invariant models for near wall turbulence. These 
two parameters will be tested in the turbulence module for the NPARC code for 
complex flow calculations. 

3.3 Wall Function Features for NPARC 

An alternative to the low-Reynolds number turbulence model approach, in which 
the integration is carried down to the wall, is to use the wall functions which provide 
boundary conditions for the solutions at the equilibrium log layer of a turbulent 
boundary layer instead of the wall. The wall function approach has the appeal that 
it reduces the near-wall grid resolution requirements considerably, and consequently 
reduces the numerical stiffness as well. However, the traditional wall functions were 
developed based on the turbulent boundary layer flow over a flat plate at zero 
pressure gradient. It does not contain the pressure gradient effect and it can not be 
used for flows with separation, a common feature for flows in propulsion systems. 
Better wall functions have been developed by Yang and Shih^^ and Shih et al.^^ 
These new wall functions take into account the effect of the pressure gradient. In 
addition, the paper by Shih et al.^^ also addresses the issue of how to implement 
wall functions to calculate flows with separation and reattachment. These new wall 
functions are going to be incorporated into the turbulence module for the NPARC 
code. We expect that the wall function feature would be available by the end of 
1995. 

3.3 Anisotropic Algebraic Reynolds Stress Model for NPARC 

All the turbulence models mentioned above are isotropic eddy viscosity models 
which assume a linear relationship between the Reynolds stresses and the mean 
strain rate. This assumption is not valid for complex flows. Anisotropic Reynolds 
stress models overcome this deficiency and will give a better performance than the 
current isotropic eddy viscosity models for flows where anisotropy is important, for 
example, the stagnation point flow and flows with massive separation. Recently, one 
such anisotropic algebraic Reynolds stresses model was developed by Shih et al.^^ 
and has shown improved predictions for a large range of complex flows including 
the confined jet flow and the swirling jet flow. This anisotropic algebraic Reynolds 
stress model is going to be implemented into the turbulence module for the NPARC 
code, and the resulting turbulence module will be used together with the NPARC 
code to calculated a large range of practical flows. To implement this model, the 
NPARC code for the mean flow solver has to be modified because in addition to 
the eddy viscosity term, the anisotropic Reynolds stresses model also brings other 
source-like terms into the mean momentum equations. The planned milestone for 
this work is the end of FY96. 
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Figure 1: Configuration for turbulent boiindary layer fiow over a fiat plate 
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Figure 2: Skin friction distribution for the turbulent boundary layer fiow 
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Figure 3: Configuration for an ejector nozzle flow 



Figure 4: Velocity distributions in the ejector nozzle flow 


ORIGINAL PAGE IS 
OF POOR QUALITY 

























Center for Modeling of Turbulence and Transition 
Research Briefs - 1995 


89 


Calculations of Turbulent Reacting 
Flows in Can Combustors 

A. T. Norris and W. W. Lion 

1. Motivation and Objective 

A newly developed probability density function (PDF) code is applied in the 
numerical simulation of swirling can combustor flows. 

Due to enviromental concern over NOi emissions, the design of future gas tur- 
bine combustors, operating at high pressures and temperatures, requires a detailed 
knowledge of the combustion processes in the different operating conditions of an 
engine. Therefore, robust and accurate computational modeling of the comubstor 
flow flelds are necessary in the design of efficient, low emission combustors. 

In modeling turbulent reactive flows, PDF methods have an advantage over the 
more traditional moment closure schemes. PDF methods treat the chemical reaction 
source term in an exact manner, while moment closure schemes model the mean 
reaction rate. For example, the commonly used laminar chemistry approximation 
is based on the assumption that the effects of turbulence on the chemistry are 
negligible. This, or similar assumptions fail for flows with high turbulence levels 
and finite rate chemistry. 

2. Hybrid PDF Method 

The hybrid PDF scheme used is a combination of moment closure and PDF 
methods. The continuity and momentum equations are solved by a moment closure 
method while the transport equation for the joint PDF of composition is used to 
solve the scalar fields. This combination has the advantage in that the PDF part 
of the code can be added to an existing moment closure method with a minimal 
amount of coding. Thus existing combustion codes can be extended to include the 
PDF model of turbulent chemical reactions without the need for extensive rewriting. 

3. Experiment 

The flow chosen to test the new code, is the axisymmetric can combustor (ASCC) 
tested by UC Irvine. [1] This combustor is cyhndrical, with an 80 mm ID and is 320 
mm long. A sketch of the combustor is shown in Fig. 1. The ASCC is operated 
at 1 atmosphere pressure. The fuel is gaseous propane, which is injected through 
an annular cone nozzle with an inner radius of 5.4 mm and an average injection 
velocity of 12.4 m/s. Air swirling is achieved by 12 curved blades. The flow rate of 
the center swriling air was kept the same as the outer dilution air. 

4. Numerical Solutions 

The mean velocity field is obtained from the solutions of Navier-Stokes equations 
using the FAST-2D code. The code uses a finite volumn scheme with cell-centered 



90 


A. T. Norris and W.W.Liou 


nodes. The pressure -velocity coupling is handled with the SIMPLEC algorithm. 
The standard high-Reynolds number k-e models were used to provide the turbu- 
lence mixing characteristics. 

The scalar PDF code is similar to that developed by Hsu et al [2] but with a several 
important numerical changes. These changes primarily involve the treatment of the 
convection and molecular mixing. Details of these changes will be discussed in a 
future paper. 

5. Initial Conditions 

The first flow measurement was reported at 0.5 cm downstream of the nozzle exit 
plane. This was found unsatisfactory for use as initial conditions for a numerical 
solution of the flow field. Therefore, the numerical initial velocity profiles for the 
streamwise, radial, and azimuthal directions were determined by the air and fuel 
mass flow rates, the angle of the swirl vanes, and the gaseus propane injection angle 
at the nozzle exit. The initial turbulent kinetic energy was set at 1% of the mean 
flow energy. The turbulent dissipation rate was then determined by assuming an 
equavilent laminar viscosity. 

Species initial conditions are simply pure fuel and air. To ignite the flame, we set 
the initial species and temperature in the combustor to correspond to a fuUy burnt 
fuel-air mixture at the experimental equivalence ratio. 

6. Thermochemistry 

Calculations were performed using the equilibrium chemistry assumption. In this 
it is assumed that the chemical time scales are much smaller than the turbulent time 
scales. In the future, a full mechanism, simplified by the Intrinsic Low-Dimensional 
Manifold (ILDM) method of Pope and Maas [3], will be used. 

7. Results 

The calculated temperature and axial velocity profiles at four different stream- 
wise locations in the combustor are compared with the measurement in Figs. 2 and 
3, respectively. For comparison, the results obtained by a moment closure model 
for the same combustor are also included. Fig. 2 shows the comparison of the axial 
mean velocity profiles. At x=l cm, the present calculation correctly predicts the 
magnitude and the radial location of the peak velocity, which roughly corresponds 
the edge of the air swirler. The reciculation zone indicated by the measurement 
near the inner edge of the swirler is not predicted well by either the present calcu- 
lation based on PDF method or the moment closure chemistry model. Therefore, 
it is possible that this may be due to the fc - e turbulence model, which is known 
to underpredict flow recirculation in highly swirling flows. At the next three axial 
stations, the PDF calculation does a better job capturing the correct local charac- 
teristics of the flow than the momentum closure procedure. 

Fig. 3 show the radial temperature distributions at the same four axial loca- 
tions as in Fig. 2. The PDF method has returned better prediction for the peak 
termperature and profiles in the region near the nozzle exit, x=l and 4 cm than 
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the momentum closure chenaistry model. The shghtly higher predicted tempera- 
ture may be due to either the insufficient turbulent mixing in the recirculation zone 
returned by the k - e model or the equilibrium chemistry. The PDF method also 
correctly predicted the centerhne temperature. At farther downstream the PDF 
model also show better temperature profile distribution. This can also be observed 
from the relatively uniform temperature coutours in this region, Fig. 4. Note that 
the experimenal flow shows higher turbulent mixing, resulting in more uniform tem- 
perature near the center region of the combustor. Again, this broad temperature 
plateau region may be due to the deficiency of the turbulence models or the equilib- 
rium chemistry. The preliminary results of finite rate chemistry calculations shows 
a lower peak temperature at x=l cm. 

8. References 

1 R. Charles, “Detailed Data Set: Velocity and temperature measurements in the 
axisymmetric can combustor (ASCC) for a parametric variation in inlet condi- 
tions” UCI-ARTR-87-6, University of California, Irvine 1988. 

2 A. T. Hsu, M. S. Raju and A. T. Norris, “Apphcation of a PDF method to 
compressible turbulent reacting flows” AIAA 94-0781, Reno NV 1994. 

3 S.B.Pope and U.Maas, “Simphfying chemical Kinetics: Trajectory-generated low- 
dimensional manifolds” FDA-93-11, Cornell University, Ithaca NY 1993. 
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Figure 1. A sketch of the UC Irvine axisymmetric can combustor 
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Figure 2. Radial distribution of the mean streamwise velocity. •, measuremen 
Micklow et. fl/[4]; , present. 





Modeling of Turbulent, Reacting Flows by PDF Methods 


93 


2500.0 

' ; 1 

' 1 


2500 .0 

: — 1 — 1 — ^ — ; — ' — 1 — ' : 

2000.0 

- 

(a) 

- 

2000.0 

(C) 

1500.0 

- r - , 

• 


- 

1500 . 0 * 


1000.0 

- 

\ 

” 

1000.0 

- 

500.0 

s 


.A 

500.0 




0.0 

■ 1 1 1 

0.0 







Figure 3. Comparison of Radial distribution of temperature. •, measurement [1]; - - 
Micklow el a/[4]; , present. 


x=1 cm 4 cm 


10 cm 


24 cm 



Figure 4. Temperature contours 







Page intentionally left blank 


Center for Modeling of Turbulence and Transition 
Research Briefs - 1995 


95 


Assessment of Turbulence 
Models in Turbomachinery Flows 


I. Motivation and Objectives 

Turbomachinery flows pose a major challenge for CFD since they are three- 
dimensional and time dependent. NASA Lewis Research Center is playing a lead 
role in developing CFD tools which can be used by industry in the analysis of their 
designs. 

Recently, a blind test case for a rotor compressor (ROTOR 37) was organized by 
the ASME/International Gas Turbine Conference held at the Hague (A. Strazisar, 

J. Wood, and K. Suder^). The objective was to judge the predictive capabilities 
of turbomachinery CFD tools. A total of 12 groups submitted their predictions 
for this case, which were then compared with the detailed experimental data. The 
results, which are given in the above reference, showed a wide scatter relative to 
the experimental data. There were several reasons for this, one of them being the 
inadequacies of the turbulence models. 

Since different codes using the same turbulence models produced different results 
in the ROTOR 37 blind test case, it is obvious that there is a need to assess 
the performance of different turbulence models from the same numerical platform. 
If different turbulence models are implemented in a single CFD code, then the 
differences in the results can be attributed only to the differences in the models. 
In this way the impact of different models on the prediction of turbomachinery 
flows can be assessed in a systematic manner. In order to answer this question, 
CMOTT has started a joint project with the turbomachinery group at the Lewis 
Research Academy. Since the Average Passage Equation code, VSTAGE, predicted 
the ROTOR 37 blind test case reasonably well^, it was selected as the numerical 
platform to carry out this task. 

The approach taken is to implement several turbulence models, of varying so- 
phistication, in the VSTAGE code in the form of a subroutine. To begin with the 
standard k-e model of Launder and Spaulding^ will be implemented. This will be 
followed by implementing an improvement to this model which has been developed 
by CMOTT. This improvement has already shown to give better results over the 
standard k — e model for a some benchmark flows. Eventually a complete CMOTT 
two equation model (Shih, et. al.^) will be implemented. 

All of the above models will be used with the wall functions. The standard wall 
functions approach in CFD calculations assumes the existence of the law of the 
wall, even though this is true for a flat plate boundary layer only. This approach 
is popular because it requires a mesh size that is computationally affordable and 
because it, apparently, provides good engineering accuracy. At a later time, and 
based on the results obtained using the above models, there may be a need to assess 
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the performance of the different wall functions on the prediction of turbomachinery 
flows. It should be pointed out that all of the turbulence model assessment wiU be 
done against the ROTOR 37 data. 

In the course of this project we may find some deficiencies in the turbulence 
models. It is also the objective of this effort to identify and remedy these. 

This project will make the following two important contributions. First an as- 
sessment of the hierarchy of turbulence models (including those that CMOTT has 
developed) in the context of turbomachinery flows. This information will be ex- 
tremely useful to the turbomachinery community. Second the models which per- 
form the best against the test cases can be incorporated in the Average Passage 
Equation codes for both the single stage and the multistage turbomachinery. This 
will lead to an improved turbomachinery CFD tool. 

2. Work in Progress 

We have finished implementing the standard two equation k — e model and its 
CMOTT improvement in the VSTAGE code. The preliminary results for the RO- 
TOR 37 test case look promising. We are currently conducting a grid sensitivity 
study. The results of this work will be presented in an upcoming report. This work 
is being carried out by A. Shabbir, J. Zhu, M. Celestina, J.J. Adamczyk, and T.-H. 
Shih. 

3. References 

^ Strazisar, A., Wood, J., and Suder, K. ASME/International Gas Turbine Confer- 
ence, The Hague, Netherlands, 1994. 

^ Celestina, M. ASME/International Gas Turbine Conference, The Hague, Nether- 
lands, 1994. 

^ Launder, B.E., and Spaulding, D.B., “The numerical computation of turbulent 
flows”, Comp. Meth. Appl. Mech. Engg., 3, 269-289, 1974. 

^ Shih, T.-H., W.W. Liou, A. Shabbir, Z. Yang and J. Zhu. “A new high Reynolds 
number k — e model” , Computers Fluids, 24, 3,227-238, 1995. 
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Appendix A 

Industry-Wide Workshop on 
Computational Turbulence Modeling 


CMOTT organized an Industry-wide Workshop on Computational Turbulence Mod- 
eling on October 6-7, 1994. The purpose of the workshop was to initiate the transfer of 
technology developed at Lewis Research Center (LeRC) to industry and to discuss the 
current status and the future needs of turbulence models in industrial CFD. The work- 
shop organizing committee and the titles of presentations are listed below. Details of the 
presentations can be found in the NASA CP 10165. 

Workshop Organizing Committee 


Industries 


C. Prakash, General Electric 
M. Sindir, Rocketdyne 
S. Syed, Pratt & Whitney 


Universities 


J.Y. Chen, University of Califomia, Berkeley 
J.L. Lumley, Cornell University 


National Aeronautics and Space Administration 

L.A. PovineUi, Chairman 
R. Mankbadi, Lewis Research Center 
D.R. Reddy, Lewis Research Center 
P. Richardson, Headquarters 
RJ. Shaw, Lewis Research Center 


Institute for Computational Mechanics in Propulsion 
T. Keith, Ohio Aerospace Institute 

A. Shabbir, Center for Modeling of Turbulence and Transition 
T.-H. Shih, Center for Modeling of Turbulence and Transition 
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Workshop Presentations 


TURBULENCE PROGRAM FOR PROPULSION SYSTEMS 

Tsan-Hsing Shih, Institute for Computational Mechanics in Propulsion and Center for ' 
Modeling of Turbulence and Transition, NASA Lewis Research Center 

TURBULENCE MODEL DEVELOPMENT AND APPLICATION AT LOCKHEED FORT 
WORTH COMPANY 

Brian R. Smith, CFD Group, Lockheed Fort Woith Company 

A SUMMARY OF COMPUTATIONAL EXPERIENCE AT GE AIRCRAFT ENGINES 
FOR COMPLEX TURBULENT FLOWS IN GAS TURBINES 
R. Zetkle and C Prakash, GE Aircraft 

THE APPL ICABILITY OF TURBULENCE MODELS TO AERODYNAMIC AND PROPULSION 
FLOWFIELDS AT McDONNELL DOUGLAS AEROSPACE 

Linda D. Krai, John A. Ladd, and Mori Mani, McDonnell Douglas Aerospace 

EXPERIENCE WITH k-e TURBULENCE MODELS FOR HEAT TRANSFER 
COMPUTATIONS IN ROTATING 

Prabhat Tekriwal, GE Corporate Research and Development - 

TURBULENCE MODELS FOR GAS TURBINE COMBUSTORS 
Andreja Brankovic, CFD Group, Pratt & Whitney 

COMBUSTION SYSTEM CFD MODELING AT GE AIRCRAFT ENGINES 

D. Burras and H. Mongia, GE Aircraft Engines, and A. Tolpadi, S. Correa, and M. Braaten, 

GE Corporate Research and Development 

RECENT PROGRESS IN THE JOINT VELOCITY-SCALAR PDF METHOD 
M.S. Anand, Allison Engine Company 

OVERVIEW OF TURBULENCE MODEL DEVELOPMENT AND APPLICATIONS AT 
ROCKETDYNE 

A.H. Hadid, E.D. Lynch, and M.M. Sindir, Rocketdyne Division, Rockwell 
International 


RECENT ADVANCES IN PDF MODELING OF TURBULENT REACTING FLOWS 
A.D. Leonard and F. Dai, CFD Research Corporation 

EXPERIENCE WITH TURBULENCE INTERACTION AND TURBULENCE-CHEMISTRY 
MODELS AT FLUENT INC. 

D. Choudhury, S.E. Kim, D.P. Tselepidakis, and M. Missaghi, Fluent Inc. 

EXPERIENCES WITH TWO-EQUATION TURBULENCE MODELS 

Ashok K. Singhal, Yong G. Lai, and Ram K. Awa, CFD Research Corporation 


PROGRESS IN SIMULATING INDUSTRIAL FLOWS USING 
CAN MORE BE ACHIEVED WITH FURTHER RESEARCH? 
Vahe Haroutunian, Fluid Dynamics International, Inc 


TWO-EQUATION MODELS: 
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TURBULENCE MODELING NEEDS OF COMMERCIAL CFD CODES: COMPLEX FLOWS 

IN THE AEROSPACE AND AUTOMOTIVE INDUSTRIES 
Bizhan A. Befrui, adapco 

TURBULENCE REQUIREMENTS OF A COMMERCIAL CFD CODE 
LP. Van Doonnaal, C.M. MueUer, and MJ. Raw, Advanced Scientific 
Computing Ltd. 

SECOND-ORDER CLOSURES FOR COMPRESSIBLE TURBULENCE 

3JL. Lumley, S. Savarese, and C.C. Volte, Mechanical and Aerospace Engineering 
Department, Cornell University 

MODELING OF TURBULENT CHEMICAL REACTION 

J.-Y. Chen, Department of Mechanical Engineering, University of California, 

Berkeley 

INTRODUCTION TO TURBULENCE SUBPROGRAM 

T.-H. Shih and J. Zhu, Institute for Computational Mechanics in Propulsion and Center for 
Modeling of Turbulence and Transition, Lewis Research Center 

DESCRIPTION OF TURBULENCE SUB-PROGRAM 

J. Zhu, Instimte for Computational Mechanics in Propulsion, NASA Lewis Research Center . 

OVERVIEW OF PROBABILITY DENSITY FUNCTION (PDF) MODELING AT LeRC 
DJL Reddy, Internal Fluid Mechanics Division, NASA Lewis Research Center 

PDF METHODS FOR TURBULENT REACTIVE FLOWS 

Andrew T. Hsu, NYMA, Inc., NASA Lewis Research Center 

A COMPOSITION JOINT PDF METHOD FOR THE MODELING OF SPRAY FLAMES 
M.S. Raju, Nyma, Inc., NASA Lewis Research Center 

IMPROVEMENTS AND NEW FEATURES IN THE PDF MODULE 
A.T. Norris, Institute for Computational Mechanics in Propulsion 
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Dr. Louis A. Povinelli 
Deputy Chief 

Internal Fluid Mechanics Division 
NASA Lewis Research Center 

CMOTT Technical Leader 

Dr. Tsan-Hsing Shih 

Senior Research Associate 

ICOMP, NASA Lewis Research Center 

Advisors 

Dr. Marvin E. Goldstein 
Chief Scientist 

NASA Lewis Research Center 
Professor John L. Lumley 
Sibley School of Mechanical and 
Aerospace Engineering 
Cornell University 


Professor Eli Reshotko 
Department of Mechanical and 
Aerospace Engineering 
Case Western Reserve University 
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Listing of current 

Members 


Names /Term 

Affiliation 

Research Areas 

Liou, William W. 
11/1990 - present 

ICOMP 

Compressible Flow 
Modeling, Weakly 
Nonlinear Wave Models 

Norris, Andrew 
4/1993 - present 

ICOMP 

PDF Turbulence Modeling, 
DNS 

Shabbir, Aamir 
5/1990 - present 

ICOMP 

Buoyancy Effects on 
Turbulence, Turbulence 
Modeling 

Shih, Tsan-Hsing 
5/1990 - present 

ICOMP 

Turbulence Modeling 

Yang, Zhigang 
7/1990 - present 

ICOMP 

Modeling of Bypass 
Transition 

Zhu, Jiang 
4/1992 - present 

ICOMP 

Application of Turbulence 
Models in Complex Flows 
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Appendix C 

List of Member’s Publications 


Barton, J. M., Rubinstein, R. and Kirtley, K. R. “Nonlinear Reynolds stress model for 
turbulent shear flows, ” AIAA Paper No. 91-0609, (Jan. 1991). 

Barton, J. M. and Rubinstein, R., “Nonlinerar algebraic Reynolds stress model for 
anisotropic turbulent flows, ” 11th U. S. National Congress of Apphed Mechanics, Tucson, 
AZ (May 1990). 

Barton, J. M. and Rubinstein, R., “Renormahzation group theory and turbulence mod- 
eling, ” 2nd International Workshop on Chaos and Turbulence, Tsukuba, Japan (Jan. 
1992). 

Barton, J. M. and Rubinstein, R., “Renormalization group analysis of turbulence- driven 
secondary flows, ” 11th Australian Fluid Mechanics Conference, Hobart, Australia (Dec. 
1992). 

Brown, S., Leibovich, S., and Yang, Z., “On the hnear instability of the Hall-Stewartson 
vortex,” Theo. and Comp. Fluid Dynamics, Vol 2, no. 1, (1990). 

Duncan, B. S., Liou, W. W. and Shih, T.-H., “A multiple- scale turbulence model for 
incompressible flow,” AIAA 93-0086 (1993). 

Duncan, B.S., Lumley, J.L., Shih, T.H. and To, W.M., ”A new model for the turbulent 
dissipation” International Conference of Fluid Mechanics and Theoretical Physics, Bejing, 
China (1992). 

Durbin, P.A. and Yang, Z., “A transport equation for eddy viscosity,” Proc. of the 1992 
Summer Program, CTR, Stanford, 1992. 

Georgiadis, N.J., Chitsomboon, T. and Zhu, J., “Modification of the two-equation tur- 
bulence model in NPARC to a Chien low Reynolds number k-e formulation,” NASA TM 
106710 (1994). 

Hsu, A.T., Raju, M.S., and Norris, A.T., Application of a PDF method to compressible 
turbulent reacting flows. In AIAA 32nd Aerospaice Sciences Meeting and Exhibit, Reno, 
Nevada, pages AIAA-94-0781, 1994. 

Hsu, A.T. and Liou, M., “Computational analysis of underexpanded jets in the hypersonic 
regime,” Journal of Propulsion and Power, 7, No. 2, (1991). 

Hsu, A.T., “Progress in the development of PDF turbulence models for combustion,” 10th 
NASP Symposium, April 23-16, 1991, Monterey, California. 

Hsu, A.T., “A study of hydrogen diffusion flames using PDF turbulence model,” AIAA 
91-1780, AIAA 22nd Fluid Dynamics Conference, June 24-26, 1991, Honolulu, Hawaii. 
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Publications List 


Hsu, A.T. and Chen, J.Y., “A continuous mixing Model for PDF simulations and its 
applications to combusting shear flows,” 8th International Symposium on Turbulence Shear 
Flows, Sept. 9-11, 1991, Munich, Germany. 

Hsu, A.T., Raju, M.S., and Norris, A.T., Application of a PDF method to compressible 
turbulent reacting flows. In AIAA 32nd Aerospace Sciences Meeting and Exhibit, Reno, 
Nevada, pages AIAA-94-0781, 1994. Huang, P.G. and Liou, W.W., “Numerical calcula- 
tions of shock-wave/boundary-layer flow interactions,” NASA TM 106694 (1994). 

Lang, N.J. and Shih, T.-H., “A critical comparison of two-equation turbulence models ” 
NASA TM 105237, (1991). 

Liou, M. S. and Steffen, Jr., C. J., “A new flux splitting scheme,” NASA TM 104404 
(1991). 

Liou, W. W. and Shih, T.-H., “Transonic flow predictions with new two-equation turbu- 
lence models,” AIAA 95-1805. 

Liou, W. W. and Shih, T.-H., “Compressible turbulent boundary layer predictions using 
a multiple-scale model,” NASA TM (in preparation) (1995). 

Liou, W. W., Shih, T.-H., and Duncan, B.D., “A multiple— scale model for compressible 
turbulent flows,” Phys. Fluids, 7 (3), 658-666 (1995). 

Liou, W. W., “Linear instability of curved free shear layer,” Phys. Fluids, 6 (2) 541-549 
(1994). 

Liou, W. W. and Shih, T.-H., “On the basic equations for sencond-order modeling of 
compressible turbulence,” NASA TM 105277 (1991). 

Liou, W. W. and Morris, P. J., “Weakly nonlinear models for turbulentmixing in a plane 
mixing layer,” Phy. Fluids, 4, 2788-2808 (1992). 

Liou, W. W. and Morris, P. J., “The eigenvalue specturm of the Rayleigh equation for a 
plane shear layer,” International Journal of Numerical Methods in Fluids, 15 1407-1415 
(1992). 

Liou, W. W., “A new energy transfer model for turbulent free shear flow,” NASA TM 
105854 (1992). 

Mansour, N N. and Shih, T.-H. and Reynolds, W. C., “The effects of rotation on initially 
anisotropic homogeneous flows,” Physics of fluid A, 3, 10, 2421-2425 (1991). 

Michelassi, V. and Shih, T.-H., “Low Reynolds number two-equation modeling of turbulent 
flows,” NASA TM 104368, (1991). . 

Michelassi, V. and Shih, T.-H., “Elliptic flow computation by low Reynolds number two- 
equation turbulence models ,” NASA TM 105376, (1991). 

Moin, P., Shih, T.-H., Driver, D. and Mansour, N., “Direct numerical simulation of a 
three-dimensional turbulent boundary layer,” Physics of Fluids A 2, 10, 1846-1853 (1990). 


Publications List 


107 


Norris, A.T. and Pope, S.B., Modeling of extinction in turbulent diffusion flames by the 
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Norris, A.T., “The application of PDF methods to piloted diffusion flames,” Ph.D. Thesis, 
Cornell University, 1993. 

Rubinstein, R. and Barton, J. M., “Renormalization group analysis of the Reynolds stress 
transport equations, nonlinear Reynolds stress models and the renormalization group,” 
Phys. Fluids A 4, 1759 1766 (1992). 
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